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Several  alternative  decouplings  of  the  electron  propagator  are 
investigated  in  this  dissertation  in  an  attempt  to  derive  more  accurate 
and  more  tractable  computational  schemes  for  extracting  the  physical 
information  contained  in  the  electron  propagator.  When  the  electron 
propagator  is  defined  as  a  single-time  Green's  function,  the  decoupling 
approximation  and  the  choice  of  reference  state  average  are  shown  to  be 
independent  approximations,  and  the  use  of  uncorrelated,  Hartree-Fock 
reference  states  is  advocated.  The  derivation  of  each  decoupling 
approximation  utilizes  the  superoperator  formalism  and  emphasizes 
elementary  algebraic  manipulations.   In  Chapter  1,  operator  product 
decouplings  are  reviewed  and  critically  discussed.   In  Chapter  2,  the 
moment  conserving  decouplings  which  consist  of  Pade'  approximants  to 
the  propagator  moment  expansion  are  investigated.   In  Chapter  3,  a 
partitioning  of  the  superoperator  Hamiltonian  is  invoked,  and  a  pertur- 
bation expansion  of  the  superoperator  resolvent  is  developed.  This 


IX 


development  leads  straightforwardly  to  the  derivation  of  the  Dyson 
equation  and  permits  an  identification  of  wave  and  reaction  superoper- 
ators.  Two  types  of  diagram  conserving  decouplings  are  then  examined, 
and  equivalences  with  the  diagrammatic  expansion  method  are  demonstrated. 
Finally  in  Chapter  4,  renormalized  decouplings  are  considered,  and  the 
two  particle-one  hole,  Tamm-Dancoff  approximation  is  specifically  derived 
and  investigated.  In  each  of  the  first  four  chapters,  the  decoupling 
approximations  are  evaluated  on  the  basis  of  computational  applications 
in  which  the  propagator  poles  are  compared  to  experimentally  determined 
ionization  energies  for  several  molecules.   In  order  to  avoid  a  possible 
bias  with  this  evaluation  criterion,  the  quality  of  the  Feynman-Dyson 
amplitudes  is  examined  in  Chapter  5  via  the  calculation  of  relative 
photoionization  intensities.  The  four  decoupling  approximations  are 
finally  summarized  as  various  approximations  to  the  wave  and  reaction 
superoperators,  and  several  extensions  of  these  investigations  are 
proposed. 


INTRODUCTION 
Since  subatomic  particles  are  beyond  the  limit  of  human  sensory 
perception,  our  knowledge  of  atomic  structure  is  based  on  the  interpre- 
tation of  measurements  with  auxiliary  probes.  The  accumulation  and 
interpretation  of  data  from  these  types  of  measurements  have  led  to  the 
conception  and  axiomatization  of  quantum  theory.  Using  the  calculus  of 
this  theory,  quantum  mechanics,  there  is  evidence  to  believe  that  it  is 
possible,  at  least  in  principle,  to  calculate  the  statistical  result  of 
any  experimental  measurement.  Unfortunately,  the  mathematical  complexity 
of  this  calculus  precludes  exact  solutions  for  all  but  a  few,  relatively 
trivial  applications;  consequently,  the  predictive  value  of  the  theory 
is  limited.  Owing  to  this  limitation,  one  aspect  of  current  theoretical 
research  involves  the  formulation  and  evaluation  of  accurate  mathematical 
approximations  which  are  relevent  to  the  interpretation  of  specific 
experiments . 

With  the  development  of  photoelectron  spectroscopy  (Turner  et  al., 
1970,  Siegbahn  et  al_. ,  1969),  photoionization  has  become  an  extremely 
useful  probe  of  atomic  and  molecular  structure  and  has  stimulated  much 
theoretical  interest  (Cederbaum  and  Domcke,  1977,  and  references  therein). 
In  the  photoionization  experiment,  light  is  shone  on  an  atomic  or  molec- 
ular sample  and  the  kinetic  energy  of  the  ionized  electrons  or  photoelec- 
trons  which  are  ejected  is  then  analyzed.  From  energy  conservation,  the 
binding  energies  of  the  photoelectrons  may  be  deduced. 


An  ab  initio,  theoretical  interpretation  of  photoelectron  spectra 
requires  the  calculation  of  ionization  energies.  These  calculations 
reflect  fundamental  assumptions  about  the  complex  nature  of  the  many- 
electron  interactions  that  occur  in  atoms  and  molecules  and  may  be  per- 
formed at  various  levels  of  sophistication.  One  conceptually  simple 
scheme  is  the  Hartree-Fock  self-consistent  field  (HF-SCF)  method  (see 
e.g.  Pilar,  1968).  In  this  method,  an  orbital  energy  is  calculated  for 
each  electron  in  an  N-electron  system  by  assuming  that  that  electron 
interacts  only  with  an  average  electron  density  formed  by  the  remaining 
N-l  electrons.  By  thus  averaging  out  the  instantaneous  electron-electron 
interactions,  the  original  N-electron  problem  is  reduced  to  N  one-electron 
problems.  The  negative  of   the  orbital  energies  obtained  in  this  calcula- 
tion can  then  be  related  to  ionization  energies  via  Koopmans'  theorem 
(Koopmans,  1933). 

Ionization  energies  obtained  at  the  Hartree-Fock  level  of  approxi- 
mation are  rarely  accurate  and  occasionally  predict  even  the  wrong  se- 
quence of  ionization.  In  order  to  obtain  more  accurate  ionization  ener- 
gies, the  electron-electron  interactions  must  be  treated  more  realisti- 
cally. In  the  Hartree-Fock  approximation,  the  N-l  electrons  of  the  ion 
state  are  assumed  to  be  "frozen"  at  the  same  energies  they  had  in  the 
ground  state.  Conceptually,  each  electron  screens  to  some  extent  the 
electrostatic  attraction  between  the  positively  charged  nuclei  and  all 
the  other  electrons  in  the  system.  As  can  be  easily  rationalized,  this 
screening  is  most  effective  for  deep-lying  or  core  electrons  which  have 
a  high  probability  density  near  the  nucleus  than  for  more  diffuse, 
valence  electrons.  Nevertheless,  if  one  electron  is  ionized,  all  the 
others  should  experience  a  stronger  nuclear  attraction  and  will  contract 


producing  a  lower  total  energy.  This  rearrangement  defines  the  relaxa- 
tion energy,  and  it  may  be  easily  incorporated  in  the  ionization  energy 
calculation.  Performing  separate  Hartree-Fock  calculations  for  both  the 
ground  and  ion  states  and  obtaining  a  total  energy  by  adding  orbital 
energies  with  corrections  for  the  overcounting  of  interelectronic  repul- 
sions, an  improved  ionization  energy  can  be  obtained  by  subtracting  total 
energies.  This  level  of  approximation  is  known  as  the  AE(SCF)  method 
(Bagus,  1965)  and  generally  yields  reliable  core  electron  ionization 
energies. 

The  remaining  discrepency  between  the  AE(SCF)  ionization  energies 
and  the  ionization  energies  obtained  from  the  exact  solution  of  a  non- 
relativistic,  many-electron  formulation  can  be  defined  as  the  correla- 
tion energy.  This  correction  arises  from  the  tendency  of  any  pair  of 
electrons  in  an  atom  or  molecule  to  correlate  their  motion  so  as  to  min- 
imize the  electron-electron  repulsion.  Electron  correlation  can  be  con- 
ceptualized as  various  virtual  scattering  events  between  bound  electrons 
in  both  the  N-  and  (N-l)-electron  systems.  The  simplest  of  these  virtual 
processes  is  a  particle-hole  excitation  in  which  one  bound  electron 
absorbs  a  virtual  photon  emitted  by  another  electron  and  is  excited  from 
its  original  Hartree-Fock  orbital  to  a  more  diffuse  orbital  of  higher 
energy.  A  hole  or  vacancy  is  simultaneously  created  and  propagates  in 
the  system.  At  some  later  time,  the  excited  electron  may  decay  re- 
emitting  the  virtual  photon  which  may  then  be  reabsorbed  by  the  first 
electron.  Adopting  the  convention  that  holes  propagate  backwards  in 
time  while  electrons  propagate  forward  in  time,  a  particle-hole  excita- 
tion in  the  (N-l)-electron  ion  state  can  be  diagrammatical ly  represented 
as 


where  the  dotted  lines  represent  the  virtual  photon  exchanges.  Although 
diagrams  of  this  type  represent  quasi-classical,  virtual  processes  which 
are  not  experimentally  observable,  they  do  provide  a  conceptual  model 
and, as  we  will  later  see  in  Chapter  3,  categorize  specific  algebraic 
expressions  that  will  be  derived  to  calculate  correlation  energy  correc- 
tions. Finally,  it  should  be  noted  that  since  the  electron  motions  will 
be  correlated  to  different  extents  in  the  N-electron  ground  state  and 
(N-l)-electron  ion  states,  it  is  not  possible  to  predict,  a  priori,  the 
effect  of  this  correction  on  the  ionization  energies  (see  Fig.  1). 

One  many-electron  formulation  which  provides  a  systematic  procedure 
for  incorporating  both  relaxation  and  correlation  corrections  into  the 
calculation  of  ionization  energies  and  which  is  the  basis  of  this  invest- 
igation, is  the  non-relativistic,  single-particle  Green's  function  or 
electron  propagator  (Linderberg  and  Ohrn,  1973).  This  method  has  the 
advantage  of  yielding  ionization  energies  directly  unlike  other  many- 
electron  formulations  which  necessitate  the  total  energy  calculation  for 
both  ground  and  ion  states  and  which  yield  the  ionization  energy  as  the 
difference.   In  the  latter  methods,  significant  loss  of  accuracy   is 
inherent  in  the  subtraction  of  two,  nearly  equal  total  energies  to  obtain 
a  much  smaller,  ionization  energy.  Care  must  also  be  taken  to  avoid 
disparate  levels  of  approximating  electron  correlation  in  the  calculation 
of  each  different  state. 
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The  electron  propagator  can  be  written  as  a  function  of  the  space 
and  spin  coordinates  of  one  electron  and  of  a  complex  energy  variable. 
For  each  bound  electron,  this  function  describes  mathematically  all  of 
the  complicated  many-electron  interactions  between  that  electron  and  the 
remaining  N-l  electrons.  This  function  also  contains  information  about 
the  interaction  of  an  additional  electron  with  the  N  bound  electrons  if 
this  electron  were  to  be  added  to  the  system  in  any  of  several  possible 
orbitals.   In  a  discrete  basis  representation  of  the  electron  propagator, 
this  information  manifests  itself  as  simple  poles  or  singularities  at 
those  values  along  the  real  energy  axis  which  correspond  to  electron  ion- 
ization energies  or  electron  affinities,  that  is,  electron  detachment  or 
attachment,  respectively.  The  residues  at  these  poles  yield  the  single- 
particle  reduced  density  matrix  from  which  the  N-electron  ground  state 
average  of  any  one-electron  operator  may  be  calculated  or  which  may  be 
related  to  transition  probabilities  for  electron  detachment  or  attachment 
provided  some  description  of  the  removed  or  added  electron  is  included. 

The  electron  propagator  can  be  calculated  in  several  ways.  The 
procedure  adopted  here  is  derived  from  the  electron  propagator  equation 
of  motion,  but  one  aspect  of  this  investigation  will  demonstrate  the 
formal  equivalence  between  this  method  and  the  diagrammatic  expansion 
technique.  The  equation  of  motion  relates  the  electron  propagator  to 
the  more  complicated  two-particle  propagator.  This  two-particle  propa- 
gator also  satisfies  an  equation  of  motion  which  relates  it  to  the  three- 
particle  propagator,  and  so  on.  This  hierarchy  of  equations  finally 
terminates  with  the  N-particle  propagator,  but  in  order  to  make  any 
practical  calculation,  this  hierarchy  must  be  approximated  at  some  lower 
level,  generally  by  expressing  an  M-particle  propagator  in  terms  of  an 


(M-l)-particle  propagator.  This  approximation  is  called  decoupling  the 
equations  of  motion  and  is  not  unique.  The  accuracy  of  calculated  ioni- 
zation energies  and  the  computational  effort  in  obtaining  them  depend 
critically  on  the  decoupling.  This  investigation  proposes  and  evaluates 
several  alternative  methods. 

In  Chapter  1,  the  electron  propagator  is  formally  defined,  and  its 
equation  of  motion  is  derived.  After  an  introduction  of  the  superoperator 
formalism,  the  electron  propagator  is  approximated  by  an  inner  projection, 
and  the  decoupling  problem  is  studied  in  terms  of  the  selection  of  an 
inner  projection  manifold.  The  remainder  of  this  chapter  discusses  the 
computational  procedure  for  solving  the  propagator  equations  and  presents 
a  critical  evaluation  of  the  operator  product  decoupling. 

Chapter  2  describes  some  general  aspects  of  the  Pade'  approximant 
method  and  its  application  to  the  calculation  of  the  electron  propagator. 
Owing  to  the  conservation  of  various  moment  matrices  in  the  propagator 
equation  of  motion  these  decouplings  are  known  as  moment  conserving 
decouplings  (Goscinski  and  Lukman,  1970).  Numerical  results  obtained 
from  the  [1,0]  and  [2,1]  Pade'  approximants  are  presented  and  discussed. 
A  partitioning  of  the  superoperator  Hamiltonian  and  a  perturbation 
expansion  of  the  superoperator  resolvent  in  the  operator  space  is  devel- 
oped in  Chapter  3.  A  superoperator  Dyson  equation  is  derived  and  wave 
and  reaction  superoperators  are  identified  in  analogy  with  ordinary 
resolvent  operator  techniques.  Truncations  of  the  wave  superoperator 
operating  on  simple  annihilation  and  creation  operators  are  shown  to 
yield  inner  projection  manifolds  that  result  in  Pade'  approximants  to 
the  self-energy.  These  Pade'  approximants  conserve  various  orders  of 
the  perturbation  expansion  for  the  self-energy  and  are  therefore 


categorized  as  diagram  conserving  decouplings.  Two  approximations  based 
on  this  decoupling  scheme  are   discussed  and  evaluated. 

Chapter  4  is  devoted  to  renormalized  decouplings.  These  decouplings 
sum  certain  types  of  self-energy  diagrams  to  all  orders.  The  two-particle, 
one-hole  Tamm-Dancoff  approximation  (2p-h  TDA),  which  sums  all  ring  and 
ladder  diagrams,  is  explicitly  derived  and  discussed  in  terms  of  the 
superoperator  formalism.  The  diagonal  2p-h  TDA  previously  proposed  by 
other  authors  (Cederbaum,  1974,  Purvis  and  Ohrn,  1974)  is  re-examined 
and  is  shown  to  neglect  certain  diagonal  contributions.  Both  approxima- 
tions are  analyzed  diagrammatical ly,  and  numerical  results  are  presented 
and  discussed. 

The  evaluation  of  each  decoupling  approximation  in  the  first  four 
chapters  is  ultimately  based  on   a  comparison  of  propagator  poles  to  ex- 
perimental ionization  energies.  Chapter  5,  on  the  other  hand,  attempts 
to  corroborate  this  evaluation  criterion  by  an  examination  of  the  quality 
of  the  Feynman-Dyson  amplitudes.  This  is  indirectly  accomplished  via 
the  calculation  of  relative  photoionization  intensities  and  their  com- 
parison with  experimental  data.  The  requisite  equations  for  the  photo- 
ionization cross-section  are  derived  in  terms  of  the  Feynman-Dyson  ampli- 
tudes, and  the  most  critical  approximations  are  discussed.  Finally, 
numerical  results  are  presented  and  are  also  discussed. 


CHAPTER  1 
OPERATOR  PRODUCT  DECOUPLINGS 

1 • 1  Definition,  Spectral  Representation,  and  Equation  of  Motion  of  the 
Electron  Propagator      

The  electron  propagator  is  most  commonly  defined  as  a  double-time 
Green's  function  which  involves  an  exact  N-electron  ground  state  average 
of  a  time-ordered  product  of  electron  field  operators,  i/j ( x , t )  and 
^+(x',t')  (Linderberg  and  Ohrn,  1973).  These  field  operators  are  gener- 
ally expressed  in  the  Heisenberg  representation, 

>Mx,t)  =  exp(iHt)  ijj(x,0)exp(-iHt)  (1.1) 

and  are   functions  of  the  combined  space-spin  and  time  coordinates  of  the 
electrons.  The  operator,  H,  in  the  exponentials  is  the  N-electron  Hamil- 
tonian.  The  field  operators  4>(x,t)  and  i|V(x',t')  have  the  property  of 
annihilating  and  creating,  respectively,  an  electron  at  the  space-spin 
and  time  coordinates,  x,t  (x',t').  Letting  k->  denote  an  exact  eigen- 
state  of  the  N-electron  Hamiltonian,  |rf_1>  and  \^+l>   denotinq  exact 
eigenstates  of  the  (N-l)-  and  (N+l)-electron,  ion  Hamiltonians,  these 
properties  are  expressed  as 

^(x,t)|^>  =  l  c.   expI-KE^  -  E^_1)t>  |vF^|-1>  (i.2) 

and 

*+(x\t')|^>  =  £  Cj  exp.{-1(Ej+1  -  E5J)t'}|^+1>  .  (1.3) 
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If  the  N-electron  Hamiltonian  is  time  independent,  it  is  easy  to  show 
that  the  double-time  Green's  function  depends  only  on  the  time  difference, 
t-t'  (Fetter  and  Walecka,  1971).  In  all  practical  calculations,  however, 
this  definition  proves  to  be  a  severe  restriction  since  the  exact  N- 
electron  ground  state  is  rarely  known  and  an  approximate  ground  state 
average  is  usually  employed.  With  an  inexact  ground  state  average,  the 
electron  propagator  will  depend  on  both  t  and  t'.  To  avoid  this  restric- 
tion, it  is  possible  to  define  the  electron  propagator  as  a  single-time 
Green's  function  by  choosing  t'  equal  to  zero  (Simons,  1976), 

<<^+(x' ,0) ;  iHx,t)»  =  -ie(t)<iHx,t)  ^'(x',o)> 

+i9(-t)<ij>+(x\0)  *(x,-t)>  •  (1.4) 

This  definition  insures  the  dependence  on  only  the  time  difference  even 
when  the  ground  state  average  is  inexact.  The  brackets,  <...>,  in 
Eq.(1.4)  represent  an  average  which  may  be  either  a  pure-state  average  or 
an  ensemble  average, 

<■■■>=  I   PK  <k\    ...    \k>   .  (1.5) 

K 

To  elucidate  the  analytic  properties  of  the  electron  propagator,  it 
is  convenient  to  derive  the  spectral  or  Lehman  representation  (Linderberg 
and  Ohrn,  1973).  This  representation  is  obtained  by  first  expanding  the 
eigenstates,  |k>,  in  terms  of  the  exact  N-electron  eigenstates,  |H'^>, 

l<>  =  *  cKi|^>  .  (1.6) 

i 

Using  resolutions  of  the  identity  in  terms  of  the  exact  (N-l)-  and  (N+l)- 
electron  eigenstates,  Eq .(1.4)  can  be  written: 
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<«J>+(x',0);  i|;(x,t)»  =  Z   Z   PK  cKi  cK-  { -  i  6  ( t )  <^^J  |^  ( x ,  t )  \^+1> 
k  i,j,k 

x  <^+V(x',0)|^N>  +  ie(-t)<^|!{)+(x1,0)|^-1><^-1|^(x,-t)|^>}  .(1.7) 

Explicitly  introducing  the  time  dependence  of  the  field  operators,  Eq .(1.7 
becomes 

«fi"(x',n);  i|)(x,t)»  =  Z   Z   P  c  .  c  . 
k   i,j,k     ]   J 

x  {-ie(t)<^V(x,0)|^+1><^+1|^(x,,0)|^N>exp[-i(E^1-E5,)t] 

+  ie(-t)<^|^t(x,,C)|l'^1><H'^1|^(x,0)|H'N>exp[i(E5,-E^1)t]}  .      (1.8) 

Since  the  calculation  of  ionization  energies  (electron  affinities)  from 
the  electron  propagator  will  be  treated  as  a  transition  between  station- 
ary states  of  the  N-electron  ground  state  and  N-l  (N+l)-electron  ion 
states,  it  is  convenient  to  Fourier  transform  Eq.(1.8)  into  an  energy 
representation, 

«^+(x');  4j(x)»e  =  f  «^+(x',0)^(x,t)»exp(iEt)dt  (1.9) 

-00 

which  yields  the  spectral  representation: 

*       r    f-il^(x)fn(x'  ) 

«^+(x');  iKx)»c  =  lim  Z   Z   P  c.c.t ^ 


"k   ui 


+    gjk(x)9kj(x')  | 

E-tEj-EJJ"1)  -  in;  ^-10) 


13 
where         f.^x)  =  <^|tHx,0)|^+1>  (1.11) 

and  g.k(x')  =  <^|^t(x-  j0)  |¥N-1>  (1-12) 

are  referred  to  as  the  Feynman-Dyson  amplitudes.  From  Eq.(l.lO)  it  is 
observed  that  the  electron  propagator  has  simple  poles  along  the  real 
energy  axis  corresponding  to  the  difference  between  exact  eigenvalues  of 
the  N-electron  Hamiltonian  and  the  (N+l)-electron  Hamil tonians.  The 
poles  of  this  function,  therefore,  have  a  physical  interpretation  as 
ionization  energies  and  electron  affinities. 

Since  atomic  and  molecular  computations  are  most  conveniently  per- 
formed in  a  Hilbert  space,  we  introduce  a  complete,  orthonormal  set  of 
one-electron  spin  orbitals,  {un-(x)}.  In  this  basis,  the  electron  field 
operators  are  represented  by  the  expansion, 

<J>(x,t)  =  I   a^tju^x)  (1.13) 

i 

which,  with  the  expansion  of  the  adjoint,  defines  the  spin  orbital  anni- 
hilation and  creation  operators,  a^t)  and  aj(0).  At  equal  times,  these 
operators  satisfy  the  usual  anti-commutation  relations, 

[ai,a+]+  =  &..  (i.i4) 

[a-,aj]+  =  [ai,aj]+  =  0  (1.15) 

In  this  discrete  representation,  the  causal  electron  propagator  can  be 
written 

«ff-(x',0);  ip(x,t)»  =  T.     u*(x')«at(0);  a.(t)»ui(x)      (1.16) 
i,j  J      J 

where 
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«at(0);  ai(t)»  =  -  ie(t)<a1  (t)a](0)>  +  ie(-t)<aj(o)ai  (-t)>  .   (1.17) 

Although  a  computational  scheme  for  obtaining  ionization  energies 
and  electron  affinities  could  now  be  established  from  the  spectral  reso- 
lution, it  is  more  convenient  to  develop  a  scheme  based  on  the  equation 
of  motion  (Linderberg  and  Ohrn,  1973), 

i  ft  «at(0);  a.(t)»  =  6(t)<a.(t)at(0)+at(0)a.  (-t)> 

+  «at(0);  [ai(t),H]_»,  (1.18) 

which  follows  directly  from  Eq.  (1.17).  The  quantity  «a;.(0);  [a.,H]  » 
is  a  two-particle  propagator,  and  the  N-electron  Hamiltonian  appearing 
in  the  commutator  has  the  following  Hilbert  space  representation: 

H=  S  h  ala     +  \  I  <rr'  |  |  ss  '>  aV.a, a         (1.19) 

r,s  rs  r  s    r,r',s,s'  r  r  s  s 

where 

hrs  =  f  u*(l)[->,V2(l)-Z  ^-]us(l)dT1  (1.20) 

'  a  la 

and 

<rr'||ss'>  =  ^   u*(l)u*1(2)r^(l-P12)us(l)us,(2)dT1dT2      (1.21) 

With  a  notation  similar  to  that  used  in  Eq.  (1.9),  the  energy  trans- 
forms of  the  various  quantities  in  Eq.  (1.18)  can  be  defined,  e.g. 


«at;  ai»E  =    «at(0);  a.(t)»  exp(iEt)dt  .  (1.22; 

-00 

Substituting  the  inverse  transforms: 

i  ft  «aj(0);a.(t)»  =  ^-  [  E«at;  ai»E  exp(-iEt)dE,      (1.23) 
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oo 

6(t)<a1(t)aT(0)+aj(0)ai(-t)>  =  ~  j  <[ai ,at]+>exp(-iEt)dE,  (1.24) 

and 

«at(0);    [a.(t),H]_»  =  ij     «at(0);    [a.  (t),H]_»Eexp(-iEt)dE, 

(1.25) 
into   Eq.  (1.18) ,  we  obtain 

00 

-M       t  t         t 

2tt  J    {E«a.;a.»E  -  <[a.,a '.]   >  -  «a   ;[a,,H]_»F}  exp(-iEt)dE  =  0.(1.26) 

_oo  «  'J  J  ' 

From  the  general  properties  of  Fourier  transforms,  it  can  be  shown 
(Morse  and  Feshbach,  1953)  that  Eq.(1.26)  implies 

E«at,a.»E  =  <[a.,at]+>  +  «at; [a1  ,H]_»£  .  (1.27) 

which  represents  the  energy  transform  of  the  equation  of  motion.  The 
iteration  of  this  equation  yields  N  coupled  equations  relating  the  single- 
particle  (electron)  propagator  to  each  of  the  higher-particle  propagators. 
Successive  substitution  of  these  more  complicated  propagators  back  into 
Eq.(1.27)  yields 

«at;ai»E  =  E"1<[ai,at]+>  +  E"2<[  [a.  ,H]_,aT]+> 

+  E"3<[[[a.,H]_,H]_,a^']+>  +  .  .  .  .  (1.28) 

1.2  The  Superoperator  Notation  and  Inner  Projection  Technique 

The  use  of  superoperators  has  antecedents  in  the  work  of  Zwanzig 
(1961)  and  Banwell  and  Primas  (1963)  in  statistical  physics  and  was  intro- 
duced into  atomic  and  molecular  propagator  theory  by  Goscinski  and  Lukman 
(1970).  As  a  notational  simplification,  the  definition  of  a  superoperator 
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Hamiltonian  and  identity,  H  and  I,  provides  a  convenient  representation 
of  the  nested  commutators  appearing  in  Eq.  (1.28).  More  formally,  this 
notation  provides  a  connection  with  the  time-independent  resolvent 
methods  introduced  into  many-body  theory  by  Hugenholtz  (1957). 

The  superoperator  Hamiltonian  and  identity  are  defined  to  operate 
on  the  spin  orbital  annihilation  and  creation  operators  through  the 
relations 


and 


Ha.  =  [ai,H]_  (1.29; 


la.  =  a.  .  (1.30; 


Powers  of  the  superoperator  Hamiltonian  are  defined  by  successive  appli- 
cation of  this  superoperator,  i.e. 

H2a.j  =  H[a.,H]_  =  [  [a.  ,H]_,H]_,  (1.31) 

and  will  always  yield  linear  combinations  of  odd  (Fermion-1 ike)  products 
of  the  simple  field  operators,  a.  and  a'.  This  set  of  all  Fermion-like 
operator  products,  {X.},  forms  a  linear  space  and  supports  a  scalar 
product  defined  by 

(Xj|X.)  =  Tr{p[X.,x]]+}  (1.32) 

where  p  is  a  normalized,  but  otherwise  arbitrary,  density  operator  corre- 
sponding to  the  N-electron  ground  state  average  of  the  electron  propagator. 

Using  the  preceding  definitions  and  notation,  Eq .(1 .28)  can  be  re- 
written as 

«at;  ai»E  =  G(E)ij  =  E"1(aj.|ai)  +  E"2(aj.|Haj) 

+  E"3(a.|H2a1)  +  .  .  .  .  (1.33) 


17 

Collecting  all  annihilation  operators  in  a  row  matrix  and  all  creation 
operators  in  a  column  and  formally  summing  the  expansion  in  Eq.  (1.33), 
the  matrix  equation  of  motion  for  the  electron  propagator  becomes: 

G(E)  =  (a|(Ei-H)_1a).  (1.34) 

The  superoperator  resolvent  in  Eq.  (1.34)  can  now  be  represented  in 
closed  form  by  a  matrix  inverse  using  the  inner  projection  technique 
(Lowdin,  1965,  Pickup  and  Goscinski,  1973).   Introducing  a  projection 
operator, 

0  =  IWIirVl  3  (1.35) 

where  0=0  and  0=0,  the  inner  projection  of  a  positive  definite, 
self-adjoint  operator,  A,  is  given  by 

A'  =  K2   6  f\h   ;       A  >  0  .  (1.36) 

Making  the  substitution 

ID  =  A"%|h)  ,  (1.37) 

the  inner  projection  of  Eq.  (1.36)  becomes  (Bazley,  1960) 

A'  =  |jl)  (ill  A"1Jh)~1(hL|  (1.38) 

and  satisfies  the  operator  inequalities  (Lowdin,  1965) 

0  <  A'  <  A  .  (1.39) 

Since  the  superoperator  resolvent  in  Eq.  (1.34)  is  an  indefinite  opera- 
tor, it  is  not  valid  to  discuss  an  inner  projection  of  the  type  in  Eq . 
(1.36).  Equation  (1.38)  however,  which  does  not  contain  A2,  is  still  an 
acceptable  definition  of  the  inner  projection  provided  A  is  nonsingular. 
Using  this  definition  for  an  indefinite  operator,  the  equality  in  Eq. 
(1.39)  will  still  hold  when  h^  is  complete,  but  the  bounding  properties 
will  now  be  lost  with  incomplete  manifolds.  Using  the  Bazley  inner 
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projection,  the  electron  propagator  has  the  following  form 

G(E)  -  U|]l)(h.|  (EI-H)h)~1(_h|a.)  (1.40) 

in  which  the  decoupling  problem  has  now  been  transformed  into  the 
problem  of  choosing  an  appropriate  inner  projection  manifold. 

1.3  The  Hartree-Fock  Propagator 

Before  proceeding  to  formulate  more  sophisticated  decoupling  schemes, 
it  is  convenient  at  this  point  to  recapitulate  the  approximations  under- 
lying all  propagator  calculations  and  to  demonstrate  the  algebraic  manip- 
ulations which  are  involved  by  examining  one  simple  decoupling  in  some 
detail.  Implicitly  assuming  the  clamped  nuclei  and  non-relativistic 
approximations,  there  are  basically  three  additional  approximations  in- 
volved in  any  scheme  for  computing  the  electron  propagator.  The  first 
is  the  truncation  of  the  complete  (infinite)  set  of  spin  orbitals, 
{u.j  (x)},  to  some  finite  subset.  This  approximation  is  also  characteris- 
tic of  the  more  conventional  wavefunction  formulations  and  has  received 
considerable  attention.  Standard  basis  sets  of  various  sizes  and  quali- 
ties  are  available  in  the  literature  (Roetti  and  Clementi,  1974,  Huzinaga, 
1965,  Dunning,  1970,  Dunning  and  Hay,  1977).  The  second  approximation 
is  the  specification  of  the  N-electron  ground  state  average  or  equiva- 
lents, a  density  operator  (Eq.  (1.32)),  in  terms  of  which  the  electron 
propagator  is  defined.  The  final  approximation  is  the  specification  of 
the  inner  projection  manifold  or  the  actual  decoupling  of  the  equations 
of  motion. 

The  simplest  approximation  to  the  inner  projection  manifold,  h,  in 
Eq.  (1.40)  is  just  the  set,  {a^.},  of  simple  field  operators.  With  this 
choice,  Eq.  (1.40)  simplifies  to 
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G(E)  =  (ll(EI-H)^)"1  (1.41) 

since 

(a|a)  =  I  (1.42) 

which  can  be  verified  by  evaluating  a  specific  matrix  element: 

(a^la.)  =  Tr{p[a.,at]+}  =  6.  .Tr{p}  =  <5..  .  (1.43) 

One  particularly  convenient  density  operator,  which  corresponds  to 
an  independent  particle,  ensemble  average,  is  the  grand  canonical  density 
operator  (Abdulnur  et  al_. ,  1972,  Linderberg  and  Ohrn,  1973), 

p  =  n  [1  -  <nk>  +  (2<nk>  -  1)  a[ak].  (1.44) 

This  density  operator  yields  the  following  results  for  various  operator 
averages: 


Tr{paras}  =  6rs<nr>  (1.45) 

Tr{pa  a, a  ,a  }  =  [6  8    ,    ,-6  ,6  ,]<n  ><n  ,>  (1.46) 

r  r  s  s    l  rs  r  s   rs  sr  J  r   r  \->--^u/ 

and  reduces  to  a  pure  state  average  when  occupation  numbers,  <n.>,  of 

zero  or  one  are  chosen. 

Considering  the  i,j-th  matrix  element  of  the  electron  propagator, 

G(E).  .  =  [E6.  .  -  (a.lHa.)]"1,  (1.47) 

the  remaining  operator  scalar  product,  (a  .  1 1  la • ) ,  can  be  evaluated  by 
first  operating  with  the  superoperator  Hamiltonian  (Eq.  1.29) 

Hai  =  [arH]_  =  l   hisas  +  h        I         <ir'  |  |  ss '  >a^,a  ,a  ,      (1.48) 
s         r '  ,  s ,  s  '  " 

then  anti -commuting  with  a.  (Eq.  1.32) 
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(a^ | Hai )  =  h..   +  I   _  <ir' | | js,>Tr{paJlas,}.  (1.49) 

Using  the  grand  canonical  density  operator  to  evaluate  the  trace 
(Eq.  1.45),  we  obtain 

(aj-|Hai )  =  h^  +  £  <is  '  |  |  js  '><n$  ,  >.  (1.50) 

The  particular  basis  of  simple  field  operators, 

ai  =  I   xikV  (1.51) 

which  satisfies  the  equation 

Hai  =  eiai  (1.52) 

diagonalizes  the  matrix,  (ajHa),  and  must  be  obtained  self-consistently. 
This  is  an  equivalent  statement  of  the  conventional  Hartree-Fock  procedure 
since  the  transformation  in  Eq.  (1.51)  also  defines  the  canonical  Hartree- 
Fock  spin  orbitals 

Ui  =  I   XikV  (I-") 

k 

The  eigenvalues,  c,  are  the  Hartree-Fock  orbital  energies. 

Substituting  these  results  into  Eq.  (1.47),  we  obtain 

G(E)..  =  (E-e^-^j  (1.54) 

for  this  simple  decoupling  scheme.  The  poles  of  this  function  occurring 
at  E  =  e.  correspond  to  the  Koopmans'  theorem  (Koopmans,  1933)  approxi- 
mation to  the  ionization  energies.  Because  of  the  analogy  between  this 
decoupling  and  the  conventional  Hartree-Fock  procedure, Eq.  (1.54)  is 
referred  to  as  the  Hartree-Fock  propagator  and  will  constitute  the  start- 
ing point  for  more  exact  approximations. 
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1.4  Operator  Product  Decoupling 

As  was  pointed  out  in  the  introduction,  the  Koopmans'  theorem  approx- 
imation to  ionization  energies  is  frequently  unreliable  and  necessitates 
the  incorporation  of  many-electron  relaxation  and  correlation  corrections. 
The  simple  decoupling  scheme  which  yielded  the  Hartree-Fock  propagator  in 
the  preceding  section  can  be  extended  to  incorporate  relaxation  and  corre- 
lation by  extending  the  inner  projection  operator  manifold,  h_,   in  Eq. 
(1.40).  One  extension  of  this  manifold,  proposed  by  Pickup  and  Goscinski 
(1973),  is  the  union  of  all  operator  subspaces  containing  different 
Fermion-like  products  of  simple  field  operators: 

h  =  {VU{4Vm}U{4akalaman}U  ■  •  ■  ■  ^-55) 

In  terms  of  the  equations  of  motion,  Eq.  (1.28),  this  type  of  decoupling 
is  equivalent  to  expressing  a  higher-particle  Green's  function  in  terms 
of  lower  ones  and  was  originally  discussed  in  the  context  of  atomic  and 
molecular  theory  by  Linderberg  and  Ohrn  (1967). 

The  formulation  of  explicit  electron  propagator  approximations  with 
this  extended  operator  product  manifold  is  simplified  by  the  use  of  the 
partitioning  technique  (LSwdin,  1962).  Having  already  derived  an  ex- 
pression for  the  Hartree-Fock  propagator  with  a  manifold  consisting  of 
simple  field  operators,  it  is  convenient  to  make  the  partition 

h  =  aUf  (1.56) 

where  a  is  the  subspace  of  simple  field  operators  and  f  represents  the 
orthogonal  complement  containing  all  higher,  Fermion-like  operator 
products.  This  partitioning  is  imposed  through  the  relations 
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(a  |  a)  =  (l\f)  =   1  ,  (1.57) 

(ill)  =  (Hi)  =  0  (1.58) 

and  leads  to  a  blocked  matrix  equation  for  the  electron  propagator: 

AaJ(EI-H)a)    -(a|Hf) 
G(E)  =  (10)      A  ~  7    )     [  "  |   .     (1.59) 

\^-(f|Ha)      (f|(EI-H)f; 

Solving  for  the  upper  left  block  of  the  inverse  matrix,  an  equation  for 
G_  (E)  is  easily  obtained 

G_1(E)  =  (aJ(EI-H)a)  -  (a  |  Hf )  (f  |  (EI-H)f  )_1(f  |  Ha)  .         (1.60) 

The  first  term  on  the  right  hand  side  of  Eq.  (1.60)  is  the  inverse 
of  the  Hartree-Fock  propagator  and  the  second  term  represents  the  relax- 
ation and  correlation  corrections  to  the  Koopmans'  theorem  ionization 
energies.  Based  on  the  resemblance  of  Eq.  (1.60)  with  the  Dyson  equation 
derived  in  quantum  electrodynamics  (Dyson,  1949),  the  Hartree-Fock  propa- 
gator can  be  identified  as  a  zeroth  order,  i.e.  uncorrected,  approxima- 
tion to  G(E)  while  the  remaining  term  is  identified  as  the  self-energy, 

G_1(E)  =  G^(E)  -  i(E)    .  (1.61) 

A  number  of  approximations  to  G(E),  based  on  different  choices  of 
the  operator  manifold  in  Eq.  (1.55),  have  been  reported  in  the  literature. 
Pickup  and  Goscinski  (1973)  chose  their  manifold  to  consist  of  single- 
and  triple-operator  products  and  replaced  the  superoperator  llamiltonian 
in  the  self-energy  by  the  Fock  superoperator  defined  by 

FX  =  [X,F]_   ;  F  =  Z   e^  .  (1.62) 
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This  approximation  was  applied  to  the  calculation  of  ionization  energies 
for  helium  and  nitrogen  by  Purvis  and  Ohrn  (1974)  and  was  later  extended 
to  include  the  full  superoperator  Hamiltonian  (Purvis  and  Ohrn,  1975a). 
Redmon  et  al_.  (1975)  have  derived  an  approximation  which  includes 
single-,  triple-,  and  quintuple-operator  products  in  h  and  have  computed 
the  ionization  energies  of  neon.  Finally,  several  approximations  have 
been  reported  using  an  inner  projection  manifold  of  single  and  triple 
products  in  conjunction  with  correlated  reference  states.  (Purvis  and 
Ohrn,  1975b,  Jtfrgensen  and  Simons,  1975). 

1.5  Method  of  Solution 

The  solution  of  Eq.  (1.60)  consists  of  finding  the  poles  and  Feynman- 
Dyson  amplitudes  of  the  electron  propagator  and  writing  a  spectral  repre- 
sentation similar  to  that  of  the  exact  propagator  in  Eq.  (1.10).  The 
procedure  for  obtaining  the  spectral  representation  from  the  Dyson  equa- 
tion has  been  discussed  by  Layzer  (1963)  and  Csanak  et  a]_.  (1971)  and 
begins  with  a  solution  to  the  energy  dependent  eigenvalue  problem: 

L(E)i(E)  =  i(E)W(E),  (1.63) 

where 

L(E)  =  (a [ Ha)  +  E(E)  ,                             (1.64) 

£(E)£+(E*)  =  1  ,  (1.65) 

and  W(E)  is  a  diagonal  eigenvalue  matrix.  Expanding  in  terms  of  the 
eigenfunctions,  £(E),  G(E)  assumes  the  form 

G(E)  =  i(E)[Ei-W(E)]"V'(E)  ,  (1.66) 
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and  the  poles  are  those  values  of  E  satisfying  the  equation, 

Ek  =  W    •  (1-67) 

The  energy  dependence  of  the  eigenvalues,  W.(E),  is  sketched  in  Fig.  2 
which  shows  that  the  poles  occur  at  the  intersections  of  these  curves 
with  a  line  of  unit  slope  passing  through  the  origin.  When  the  inner 
projection  manifold  is  energy  independent,  the  slopes  of  the  W,(E)  curves 
are  always  negative  since 

3E  2(E)  =  -(a|Hf)(f|(EI-H)f)-2(f|Ha)  ,  (1.68) 

and  the  number  of  propagator  poles  between  any  pair  of  self-energy  poles 
is  equal  to  the  number  of  basis  functions  in  that  symmetry  (Purvis  and 
Ohrn,  1974). 

From  the  spectral  representation,  it  was  noted  that  the  exact  propa- 
gator has  only  simple  poles,  and  it  is  easily  shown  that  the  residues  at 
the  poles  are  precisely  the  Feynman-Dyson  amplitudes.  Assuming  that  the 
approximate  propagator  in  Eq.  (1.60)  also  has  only  simple  poles,  the 
residues  can  be  obtained  from  elementary  residue  calculus  as 

lirn  (E-Ek)Gij(E)  =  ^(E^^)  (1.69) 

where 

Fk1=  t1  -dTVE^E,  •  (1-70) 

k 

According  to  the  Mi ttag-Leffler  theorem  (Mi ttag-Leffler,  1880),  the  elec- 
tron propagator  can  now  be  written  as 


Figure  2.  A  sketch  of  the  energy  dependence  of  the  function 
WijE)  between  self-energy  poles  (indicated  by 
vertical  dashed  lines).  Propagator  poles  occur 
at  the  intersections  of  these  curves  with  a  line 
of  unit  slope. 
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r  *  (e  H*  (e  ) 

Gu(E)  =  I      (E-ilk)J  (1-71: 


which  has  the  form  of  the  spectral  representation. 

1.6  Analysis  and  Limitations  of  the  Operator  Product  Decoupling 

In  order  to  analyze  approximate  electron  propagator  expressions  as 
obtained  in  Eq.  (1.60),  it  is  necessary  to  compare  the  propagator  poles 
with  the  exact  energy  differences  between  the  corresponding  M-  and  N-l 
(N+l)-electron  states.  This  analysis  can  be  performed,  in  principle, 
in  one  of  two  ways.  Since  a  full  configuration  interaction  (CI)  calcu- 
lation will  yield  the  exact  total  energy  to  any  finite  dimensional,  model 
problem,  the  total  energies  of  the  N-  and  N-l  (N+l)-electron  state  could 
be  calculated  and  then  subtracted  to  yield  the  ionization  energy  (electron 
affinity).  A  full  CI,  however,  is  not  practical  except  for  systems  well 
described  by  small  basis  sets  (<10)  because  the  number  of  configurations 
in  the  CI  expansion,  given  by  Weyl's  formula  (Shavitt,  1977),  increases 
roughly  as  N"  (2Me/N)  where  M  is  the  number  of  electrons,  M  is  the  size 
of  the  spin  orbital  basis,  and  e  is  the  constant  2.718.  Furthermore,  the 
CI  solution,  which  is  expressed  as  a  determinental  expansion,  is  not 
readily  amenable  to  detailed  analysis.  On  the  other  hand,  a  perturba- 
tion expansion  of  the  N-  and  N-l  (N+l)-electron  total  energies  also 
yields  the  exact  solution  (in  a  non-zero  region  of  convergence),  but  in 
addition,  allows  an  order  by  order  analysis  of  the  total  energy  contribu- 
tions in  terms  of  orbital  energies  and  two-electron  integrals. 

Using  Rayleigh-Schrodinger  perturbation  theory  (RSPT)  to  represent 
the  N-  and  (N-l)-electron  states  through  second  order,  Pickup  and 
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Goscinski  (1973)  derived  a  second-order  electron  propagator  expression 
for  the  energy  difference.  Since  a  Hartree-Fock  self-consistent  field 
solution  was  assumed  as  the  unperturbed  problem,  the  orbital  basis  sets 
of  the  N-  and  (N-l)-electron  states  are  different.  Before  the  total  en- 
ergy expressions  may  be  subtracted,  therefore,  the  orbitals  of  the  (N-l)- 
electron  problem  must  be  expanded  in  terms  of  the  N-electron  orbitals. 
This  procedure  has  been  extended  by  Born  et  aj_.  (1978)  through  third 
order,  and  the  resulting  third-order  self-energy  is  listed  in  Appendix  1. 
Each  third-order  term  in  Appendix  1  is  characterized  by  a  diagram  and 
may  be  alternatively  obtained  using  diagrammatic  techniques,  but  for 
pedagogic  reasons  this  discussion  will  be  deferred  until  Chapter  3. 

Although  the  results  of  Purvis  and  Ohrn  (1974,  1975a)  and  Redmon 
et  a]_.  (1975)  represent  significant  improvements  to  the  Koopmans'  theorem 
and  AE(SCF)  approximations  for  the  ionization  energies  they  computed,  the 
operator  product  decoupling  was  demonstrated  to  have  certain  computation- 
al and  formal  limitations  or  ambiguities.  Computationally,  the  most 
severe  limitation  is  the  dimension  of  the  inverse  matrix  in  the  matrix 
product  of  Eq.  (1.42).  The  dimension  of  this  matrix  increases  rapidly 
with  the  size  of  the  spin  orbital  basis  and  prohibits  an  exact  inversion 
for  all  but  the  smallest  basis  sets.   In  contrast  to  the  CI  matrix  where 
only  the  lowest  few  eigenvalues  of  each  symmetry  are  usually  computed, 
the  inversion  of  this  matrix  requires  all  the  eigenvalues  and  eigen- 
vectors. As  an  illustration  of  the  size  problem,  when  £  consists  of  only 
triple  products,  the  dimension  is  roughly  proportional  to  NM(M-N)  where 
N  is  the  number  of  electrons  and  M  is  the  size  of  the  spin  orbital  basis. 
This  limitation  necessitates  the  approximation  of  the  inverse  matrix  in 
diagonal  or  near-diagonal  form  (Purvis  and  Ohrn,  1974). 
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The  formal  limitations  of  the  operator  product  manifolds  involve 
certain  ambiguities  when  extending  the  decoupling  approximation  and  the 
difficulty  in  performing  an  order  analysis.  Although  Manne  (1977)  has 
proven  that  the  set  of  all  Fermion-like  operator  products  forms  a  com- 
plete set,  the  extension  to  higher  operator  products  is  not  consistent 
with  an  order  by  order  extension  of  the  perturbation  analysis.  In  fact, 
Redmon  et  al_.  (1975)  have  suggested  that  perhaps  some  quintuple  products 
should  be  preferentially  included  before  all  triple  products.  This,  they 
felt,  was  particularly  important  in  describing  shake-up  processes,  i.e. 
ionization  plus  a  simultaneous  excitation  of  the  (N-l)-electron  ion. 
This  observation  has  recently  been  confirmed  by  Herman  et  aJL  (1978)  in 
calculations  employing  the  closely  related  equation  of  motion  (EOM) 
method  (Rowe,  1968,  Simons  and  Smith,  1973). 

As  mentioned  in  Section  1.3,  correlation  corrections  may  be  included 
in  either  the  density  operator  or  the  inner  projection  manifold.  This 
dichotomy  leads  to  another  ambiguity:  should  larger  operator  products  be 
chosen  to  extend  the  decoupling  or  a  more  highly  correlated  reference 
state?  It  has  been  shown  by  J0rgensen  and  Simons  (1975)  that  in  order 
to  obtain  a  decoupling  approximation  correct  through  third  order,  an 
inner  projection  manifold  consisting   of  single  and  triple  products  must 
be  chosen  as  well  as  a  reference  state  which  includes  all  single  and 
double  excitations.  Unfortunately,  this  combination  of  approximations 
makes  the  order  analysis  unnecessarily  complicated,  as  we  will  later 
show  in  Chapter  3. 

Finally,  the  Hermiticity  problem  should  be  mentioned.  With  a  den- 
sity operator  that  commutes  with  the  Hamil tonian,  the  Hermiticity  of  the 
superoperator  Hamil tonian  can  be  expressed  as 
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(Xi |HXj)  =  (HX-IXj)  =  (Xj |HX1)  .  (1.72) 

With  density  operators  that  do  not  commute  with  the  Hamiltonian,  however, 

Eq.  (1.72)  is  generally  not  satisfied  and  leads  to  both  computational 

and  formal  complications.  This  problem  has  been  studied  by  Nehrkorn 

et  a]_.  (1976)  who  observed  computationally  that  the  non-Hermitian  terms 

which  arise  when  the  density  operator  is  correlated  to  first  and  second 

order  in  RSPT  are  cancelled  when  the  reference  state  was  improved  to 

second  and  third  order  respectively.  A  general  proof  was  given  by 

* 
Linderberg  which  states  that  the  Hermiticity  error  is  of  order  n+1  when 

the  reference  state  is  correlated  through  order  n. 

In  the  following  chapters,  alternative  decoupling  procedures  will 

be  proposed  and  investigated  with  the  intention  of  remedying  the  various 

limitations  inherent  in  the  operator  product  decoupling  as  discussed 

here,  yet  which  retain  a  quantitative  description  of  ionization  processes, 


'See   Ref.  (14)  in  Nehrkorn  et  al.  (1976). 


CHAPTER  2 
MOMENT  CONSERVING  DECOUPLINGS 

2.1  Pade'  Approximants  and  the  Extended  Series  of  Stieltjes 

The  evaluation  of  special  functions  assumes  a  central  role  in 
applied  mathematics.  A  large  number  of  these  functions,  from  the  simple 
trigonometric  and  exponential  functions  to  the  more  complex,  hypergeomet- 
ric  functions  and  Green's  functions,  have  power  series  expansions.  Their 
evaluation,  therefore,  consists  of  summing  the  corresponding  series  expan- 
sion. When  the  series  is  slowly  convergent  or  when  only  a  limited  number 
of  expansion  coefficients  are  known  (as  e.g.  through  perturbation  theory), 
it  may  not  be  practical,  or  even  possible,  to  evaluate  the  series  term  by 
term  until  a  desired  accuracy  has  been  achieved.  In  these  cases,  optimal 
approximations  based  on  a  limited  number  of  expansion  coefficients  are 
sought.  This  general  problem  was  first  studied  by  Tchebychev  (1874)  and 
Stieltjes  (1884)  for  the  series  which  bear  their  namesand  is  referred  to 
as  the  problem  of  moments.   (For  more  recent  reviews  of  this  problem  see 
e.g.  Wall,  1948,  Shohat  and  Tamarkin,  1963,  or  Vorobyev,  1965).  A  general 
solution  of  this  problem  was  given  by  Pade'  (1892)  and  is  known  as  the 
Pade'  approximant  method  (Baker,  1975). 

Given  a  function,  f(z),  (z  complex)  which  admits  the  formal,  but 
not  necessarily  convergent,  power  series  expansion 


Hz)  =   ?  a  zk  ,  (2.1! 

k=0  K 
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the  [N,M]  Pade'  approximant  is  defined  as  a  rational  fraction  of  the 
form  P(z)/Q(z)  where  P(z)  is  a  polynomial  of  degree  M  and  Q(z)  is  a 
polynomial  of  degree  N.  The  coefficients  of  these  polynomials  are 
uniquely  determined  by  equating  like  powers  of  z  in  the  equation 

f(z)Q(z)-P(z)=0  (through  order  zN+M)  (2.2) 

with  the  auxiliary  condition  Q(0)=1.  The  expansion  of  P(z)/Q(z),  there- 
fore, coincides  with  Eq.  (2.1)  through  the  (N+M)-th  power  of  z  and  pro- 
vides an  approximation  to  the  remaining  terms. 

The  term  by  term  convergence  of  Eq.  (2.1)  is  replaced  by  the  con- 
vergence of  sequences  of  approximants  (such  as  [N,N],  N=l,  2,  3,  ...  ) 
in  the  Pade'  approximant  method,  and  although  general  convergence 
theorems  are   difficult  to  prove  for  arbitrary  series,  there  exist  several 
extensive  special  cases  for  which  convergence  has  been  proven.  For  these 
series,  the  Pade'  approximant  can  often  be  shown  to  extend  the  natural 
region  of  convergence  (Baker,  1975)  and  may  be  viewed  as  a  method  of 
approximate  analytic  continuation.  A  sequence  of  Pade'  approximants, 
therefore,  may  converge  rapidly  when  the  original  series  expansion  con- 
verges slowly  or  not  at  all. 

Two  series  which  have  been  extensively  studied  in  the  problem  of 
moments  and  for  which  sequences  of  Pade'  approximants  have  been  proven 
to  converge  are  the  series  of  Stieltjes  (Stieltjes,  1894)  and  the  ex- 
tended series  of  Stieltjes  (Hamburger,  1920,  1921a,  1921b,  also  known  as 
the  Hamburger  moment  problem).  A  series  is  of  the  Stieltjes  type  if  and 
only  if  the  coefficients,  a,  in  Eq.  (2.1),  can  be  identified  as  moments 
of  a  distribution 
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ak  =  I  x  #(x) 


!2.3: 


where  ijj(x)  is  a  bounded,  non-decreasing  function  with  infinitely  many 
points  of  increase  in  the  interval  [o,»).  The  extended  series  of 
Stieltjes  is  defined  similarly  for  the  extended  interval  (-~,»). 

The  extended  series  of  Stieltjes  has  particular  significance  owing 
to  its  intimate  relationship  with  resolvents  of  Hermitian  operators. 
For  any  operator,  A,  we  can  define  the  operator  function 


R(zA)  =  (i-za; 


1 


!2.4) 


which  is  trivially  related  to  the  resolvent  of  A.  When  A  is  Hermitian, 
the  spectral  theorem  (Riesz  and  Sz.-Nagy,  1955)  insures  a  unique  integral 
representation  of  R(zA)  having  the  form 


R(zA) 


dE(A) 
1-zA 


[2.5] 


The  operator  E(X)  is  called  an  orthogonal  resolution  of  the  identity, 
and  when  A  has  only  a  discrete  spectrum,  it  can  be  written 


E(A)  =  E  e(A-ak)|<j>k>«J>k| 


(2.6) 


where  an  and  <j>  are   the  eigenvalues  and  eigenfunctions  of  A.  For  any 
vector  f  in  the  domain  of  An  for  all  n,  the  function 


<f|R(zA)f>  =  R  (z) 


dEf(x; 

~T^zT 


!2.7) 


represents  either  an  extended  series  of  Stieltjes  or   a  rational  fraction 
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depending  on  whether  Ef(A)  has  an  infinite  or  finite  number  of  points  of 
increase  (Masson,  1970). 

In  view  of  possible  applications  of  the  Pade'  approximant  method  to 
the  superoperator  resolvent,  we  state  two  theorems  regarding  the  extended 
series  of  Stieltjes  and  discuss  some  properties  of  two  particular  se- 
quences of  Pade'  approximants  to  these  series. 

Theorem  1:   (Wall,  1948,  theorem  86.1)  A  necessary  and  sufficient 
condition  for  f(z)  to  be  an  extended  series  of  Stieltjes  is 


det 


aQ  ax 
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.  a 
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n+1 


a  a  . i  .  .  .  a0 

n  n+1       2 


-/ 
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;  n=0,  1,  2,  . 


!2.8] 


Theorem  2:   (Masson,  1970,  theorem  4)  If  f(z)  is  an  extended  series 
of  Stieltjes  and  the  associated  moment  problem  is  determinant*, 
then,  for  fixed  j=0,+l,+2,  .  .  .  +m,  the  sequence  [N,N+2j+l] 
of  Pade'  approximants  converges  to  f(z)  for  Im  {z)f   0.  The 
convergence  is  uniform,  i.e. 


lim  ||[N,N+2j+l]-f(z)||  =  0 


:2.9! 


with  respect  to  z  in  any  compact  region  in  the  upper  or  lower 
half-z  plane. 

In  addition  to  being  uniformly  convergent,  sequences  of  [N,N]  and 

[N,N-1]  Pade'  approximants  to  extended  series  of  Stieltjes  have  two 

other  features  which  make  them  particularly  attractive  for  computational 

applications.  First,  these  approximants  are  closely  related  to 


*The  moment  problem  is  said  to  be  determinant  if  there  is  a  unique, 
bounded,  non-decreasing  function  ip ( x )  satisfying  the  moment  conditions 
in  Eq.  (2.3)  and  the  supplementary  conditions  i|j(-°°)=0  and 

ij»(x)  =  lim  k{\p(x+e)   +  ip(x-e)}  . 

e+0 
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variational  methods  (Nuttall,  1970,  1973).  When  the  operator  R(zA)  (de- 
fined in  Eq.  (2.5))  is  positive  definite,  the  [N,N]  and  [N.N-1]  approxi- 
mants  provide  the  following  bounds  to  Rf(z)  (Goscinski  and  Brandas,  1971): 
[N,N]  >  Rf(z)  >  [N.N-1]  .  (2.10) 

For  resolvent  operators  such  as  the  superoperator  resolvent  which  are 
indefinite,  bounding  properties  are   more  difficult  to  establish.  Vorobyev 
(1965)  has  shown,  however,  that  the  inverse  poles  of  the  [N,N-1]  approx- 
imant  to  R^(z)  are  equivalent  to  the  eigenvalues  obtained  from  the  usual 
Rayleigh-Ritz  variational  problem 

where  ip  =  cQf  +  c,Af  +  .  .  .  cN_,A  "  f,  and  the  coefficients  {c}  are 
variationally  determined.  In  this  sense,  the  poles  of  [N,N-1]  to  Rf(z) 
are  variationally  optimum,  but  they  have  no  definite  bounding  properties. 
The  second  attractive  feature  of  the  [N,N]  and  [N,N-1]  approximants 
is  the  ease  with  which  they  may  be  computed.  Rather  than  solving  Eq. 
(2.2)  to  obtain  the  coefficients  of  the  polynomials  P(z)  and  Q(z),  these 
approximants  may  be  expressed  directly  in  terms  of  the  series  coefficients 
{a, }  using  matrix  formulae  derived  by  Nuttall  (1967)  and  Goscinski  and 
Brandas  (1971).  For  the  [N,N-1]  approximant,  we  have 

[N.N-1]  =  a^-z/y"1^  ,  (2.12) 

where,  in  general,  a.  is  a  column  vector  with  the  elements  a.,  a.,,, 
— i  i   i+l 

.  .  .  a.+N  ,,  and  A.  is  an  N  x  N  square  matrix  with  the  columns  a_. , 
a_.  +  1,  .  .  .  a-+N_..  Similarly  for  the  [N,N]  approximant  we  can   write 

[N,N]  =  aQ  +  zajt^-z^]"1^  .  (2.13) 
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2.2  Moment  Conserving  Decoupling 

Expanding  the  superoperator  resolvent  in  Eq.  (1.34)  and  multiplying 
both  sides  of  the  equation  by  E,  the  electron  propagator  can  be  expressed 
as  the  moment  expansion 


EG(E)  =  Z  E"k(a|Hka 
k=0 


(2.14) 


Before  the  Pade'  approximant  method  may  be  applied  to  this  equation,  how- 
ever, the  conventional  definition  of  the  Pade'  approximant  must  be  gener- 
alized to  matrix  Pade'  approximants  (Baker,  1975).  This  generalization 
is  achieved  by  replacing  the  moment  coefficients  by  the  corresponding 
moment  matrices  and  noting  that  these  matrices  do  not  commute  when  per- 
forming subsequent  algebraic  manipulations.  Using  Eq.  (2.12)  to  repre- 
sent the  [N,N-1]  approximant  to  EG(E),  we  obtain 


EG(E)  =  ^(Mq-E"1^)"1^ 


(2.15) 


or  multiplying  each  side  of  this  equation  by  E  ,  Eq.  (2.15)  becomes 


G(E)  =  ^(EMq-Mj)"1^ 


(2.16! 


where  m~  is  now  a  column  matrix  with  block  elements 


% 


M 


^N-l 


;  £.  =  (al^a) 


If  c.  has  the  dimensions  M  x  M,  M.  is  an 

columns  m. ,  m. , , ,  .  .  .  m. ...  , . 
-l  -l+l       -i+N-1 


(2.17 


x  NM  square  matrix  with 
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There  is  a  close  relationship  between  Eq.  (2.16)  and  the  inner  pro- 
jection of  the  superoperator  resolvent 

G(E)  =  (a|h)(h|(Ei-H)h)_1(h|a)  .  (2.18) 

Goscinski  and  Lukman  (1970)  have  shown  that  if  the  inner  projection  man- 
ifold is  chosen  to  consist  of 

h  =  {ak}U{Hak}U  .  .  .  U{HN_1ak}  ,  (2.19) 

the  inner  projection  and  the  [N.N-1]  Pade'  approximant  are  equal.  Since, 
in  general,  the  [N,M]  Pade'  approximant  conserves  the  first  N+M+l  moments 
in  the  moment  expansion,  this  choice  of  inner  projection  manifolds  for 
the  superoperator  resolvent  is  called  a  moment  conserving  decoupling  of 
the  equation  of  motion. 

An  examination  of  the  [N,N-1]  approximant  to  the  electron  propagator 
shows  that  its  poles  are  given  by  the  eigenvalues  of 

M.c  =  MQcd  .  (2.20) 

The  matrix  NL  corresponds  to  a  metric  matrix  and  by  virtue  of  the  opera- 
tor scalar  product,  is  always  positive  definite 

Mq  =  (h|h)  =  Tr{p[hh++h+h]}>0  .  (2.21) 

The  determinants  of  the  metric  matrices  corresponding  to  various  trunca- 
tions of  the  moment  conserving  inner  projection  manifold,  i.e. 

det  (h^)  >  0        h^  =  {a}  (2.22) 

det  (hjlhj)  >  0        hj  =  {a_}U{Ha_}      ,  (2.23) 
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provide  the  necessary  and  sufficient  conditions  of  Theorem  1  to  prove 
that  the  electron  propagator  is  an  extended  series  of  Stieltjes.  Con- 
sequently, the  sequence  of  [N,N-1]  Pade'  approximants  to  the  electron 
propagator  should  be  uniformly  convergent  and  should  have  variational ly 
optimum  properties. 

The  spectral  representation  of  the  electron  propagator  (Eq.  (1.10)) 
consists  of  two  summations,  one  which  has  poles  in  the  lower  half  of  the 
complex  E-plane  corresponding  to  ionization  energies  and  one  which  has 
poles  in  the  upper  half  plane  corresponding  to  electron  affinities. 
Based  on  the  physical  argument  that  the  removal  of  an  electron  from  a 
stable  atomic  or  molecular  system  always  requires  energy,  we  might  sus- 
pect a  separation  of  the  superoperator  resolvent  which  yields  a  nega- 
tive definite  operator  for  these  processes.  If  this  were  possible,  the 
poles  of  the  [N , N-l ]  approximants  would  then  be  upper  bounds  to  the  exact 
ionization  energies  obtainable  with  a  given  basis.  This  separation 
has  not  been  explicitly  demonstrated  but  an  overwhelming  amount  of 
numerical  data  seemsto  substantiate  this  conjecture.   In  particular,  the 
[1,0]  approximant,  which  is  easily  verified  to  be  the  Hartree-Fock  propa- 
gator, generally  yields  poles  larger  in  absolute  value  than  experimental 
ionization  energies.  One  possible  exception  to  this  rule  may  be  the 
near  Hartree-Fock  limit  calculation  of  Cade  et  ah  (1966)  on   diatomic 
nitrogen.   In  this  calculation,  the  magnitude  of  the  Itt  orbital  energy 

is  slightly  (0.2  eV)  smaller  than  the  experimental  Itt  ionization  energy. 

2 
If  on  the  other  hand,  the  X  tt  ion  state  was  fortuitously  better  described 

u  J 

than  the  ground  state  with  their  extended  basis,  this  result  may  still  be 
an  upper  bound  to  the  exact  ionization  energy  i_n  that  basis. 
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Relaxation  and  correlation  corrections  are  incorporated  in  any 
[N,N-1]  Pade'  approximant  beyond  the  [1,0]  or  Hartree-Fock  approximant. 
In  particular,  we  have  studied  the  [2,1]  approximant  in  some  detail. 
This  approximant  corresponds  to  the  truncation 

h  =  {a_}U{Ha_}  (2.24) 

of  the  inner  projection  manifold  and  conserves  the  first  four  moment 
matrices.  The  operators  {f .  If^Ha.},  which  were  evaluated  in  Eq.  (1.48), 
consist  of  a  sum  over  all  triple  products  of  simple  field  operators  with 
each  operator  product  in  the  sum  weighted  by  an  antisymmetrized,  two- 
electron  integral.  These  linear  combinations  provide  a  significant  re- 
duction in  the  subspace  of  triple  operators  thus  overcoming  one  major 
limitation  of  the  operator  product  decoupling. 

Another  type  of  moment  conserving  decoupling  of  the  electron  propa- 
gator equations  of  motion  has  been  analyzed  by  Babu  and  Ratner  (1972). 
This  decoupling  is  achieved  by  truncating  the  moment  expansion  after  the 
m-th  moment  and  replacing  the  m-th  moment  matrix  with 

c  =  c  G(E)   .  (2  25) 

Solving  the  truncated  moment  expansion  for  G(E)  yields 

(Em+1I-c  )G(E)  =V  Em"kc,  (2.26) 

~™      k=0    -* 


G(E)  =  (E'^I-c)-1  'V  Em-kc,  (2.27) 

'"    k=0    ~* 

These  rational  approximants  formally  conserve  m  moments  but  are   not  of 

the  Pade1  type  (as  Babu  and  Ratner  incorrectly  identify  them)  since  the 
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auxiliary  equation,  Q(0)  =  1,  is  not  satisfied,  i.e. 

Q(0)  =  -c^  f  1     .  (2.28) 

The  auxiliary  equation  guarantees  the  uniqueness  of  the  Pade1  approx- 
imants:  only  one  [N,N-1]  Pade'  approximant  will  conserve  exactly  m 
moments.  The  nonuniqueness  of  Babu  and  Ratner's  decoupling  scheme  is 
easily  demonstrated  by  replacing  Eq.  (2.25)  with 


m   -k     m    k 

Z  E  V  =  I     E  Kc,  G(E)    0  <  n  <  m  .  (2.29) 

k=n    "^   k=n    "* 


Solving  for  G(E) , 


m  n-1 

[El  -  E  E  Kc, ]  G(E)  =  E  E'c,  (2.30] 

k=n    "*        1=0    -l 


G(E)  =  [EmI-  I     Em\]-1  V  Em-\,    ,  (2.31) 

k=n    "*    1=0      ' 


we  obtain  m  rational  approximants  (n-1,  .  .  .  m)  which  formally  conserve 
m  moments.  Because  these  approximants  are  not  uniquely  defined,  we  will 
only  consider  Pade'  approximants  in  this  chapter. 

2 . 3  Method  of  Solution 

The  first  step  in  obtaining  the  spectral  representation  of  the  elec- 
tron propagator  with  the  moment  conserving  decoupling  is  the  evaluation 
of  the  necessary  moment  matrices.  The  first  four  moment  matrices  which 
are  necessary  to  construct  the  [2,1]  approximant  have  been  evaluated  by 
Redmon  (1975)  using  the  grand  canonical  density  operator  (Eq.  (1.44)). 


An  independent  check  of  these  derivations,  however,  revealed  an  error 
in  the  matrix  elements  of  (ajH  a)  (Redmon,  1975,  Eq.  (11.30)).  The 
correct  result  has  subsequently  been  verified  by  Redmon  (1977)  and 
appears  in  Appendix  2. 

Once  the  moment  matrices  have  been  evaluated,  the  matrices  NL  and 
M,  are  constructed  and  the  corresponding  eigenvalue  problem,  Eq.  (2.20), 
must  be  solved.  In  general,  the  dimension  of  the  eigenvalue  problem  in- 
creases linearly  with  the  size  of  the  inner  projection  manifold,  i.e.  the 
[N,N-1]  approximant  presents  an  eigenvalue  problem  of  dimension  NM  where 
M  is  the  size  of  the  spin  orbital  basis.  For  the  [2,1]  approximant, 
therefore,  the  dimension  of  this  problem  is  only  twice  the  size  of  the 
spin  orbital  basis.  This  means  that  for  even  rather  large  basis  sets, 
standard  matrix  eigenvalues  techniques  may  be  employed  to  solve  this 
problem  in  nonpartitioned  form.  As  a  consequence,  all  the  poles  and 
the  spectral  density  of  the  electron  propagator  are  obtained  from  a 
single  matrix  diagonalization  thus  avoiding  the  energy-dependent  pole 
search. 

Denoting  the  eigenvectors  by  c,  where 

cc+  =  Ml1  ,  (2.32) 

and  the  eigenvalues  by  the  diagonal  matrix  d_,  the  spectral  representation 
of  the  electron  propagator  can  be  derived, 

G(E)  =  mj  cc'1[EM0-M1]"1M0  cc"1^"1^  (2.33) 

=  mj  cIEl-dfV1^1!^  .  (2.34) 

Defining  the  matrix 

2L  =  Eq  £  '  (2.35) 


42 

which  is  rectangular  with  the  dimensions  M  x  NM,  Eq.  (2.34)  becomes 

6(E)  =  x(EI-d)_1x+  .  (2.36) 

This  equation  conserves  the  first  2N  moment  matrices  of  the  moment  expan- 
sion which  implies,  in  particular, 

xx  =  (aja)  =  I  (2.37) 

from  the  conservation  of  the  first  moment. 

The  complete  solution  of  the  electron  propagator  which  is  conven- 
iently obtained  with  this  decoupling  can  be  used  to  determine  a  self- 
consistent,  single-particle  reduced  density  matrix  (1-matrix).  The  i, 
j-th  element  of  the  1-matrix  can  be  computed  from  the  contour  integral 
(Linderberg  and  Ohrn,  1973) 

<aT3j>  =  (2tt.)_1  |  G(E)ijdE  .  (2.38) 

The  contour,  c,  runs  from  -«  to  °°  along  the  real  axis  and  encloses  only 

poles  of  G(E)..  that  lie  below  the  chemical  potential  (u)  when  finally 

closed  in  the  upper  half  of  the  complex  E-plane.  The  integral  is  then 

evaluated  using  the  Cauchy  residue  theorem  (Morse  and  Feshbach,  1953) 

<a'}a.>  =  Z       lim  (E-dk)G(E)..  (2.39) 

J    k<u  E-xi,,         1J 


=  I    xx   .  (2.4o; 


Owing  to  the  orthonormal ity  of  the  spectral  density  elements  (Eq.  (2.37)) 

*      i,2 
x.,  x ..  =  6 .  .  x.,    ,  (2  41) 

ik  jk    ij  '  lk1  v^--+i; 
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it  follows  that  the  1-matrix  is  diagonal  with  occupation  numbers  deter- 
mined by 


<n,>  -  £   x. 


2 


k<u 


ikl    •  (2.42) 


Using  pure  state  occupation  numbers  of  zero  and  one  in  the  grand  canoni- 
cal density  operator  for  the  initial  computation  of  G(E),  then  occupation 
numbers  determined  from  Eq.  (2.42)  on  subsequent  computations,  a  self- 
consistent  set  of  occupation  numbers  can  be  sought. 

2.4  Computational  Considerations  and  Applications 

The  most  time  consuming  step  in  the  construction  of  the  [2,1]  Pade' 
approximant  to  the  electron  propagator  is  the  construction  of  the  moment 
matrices.  The  fourth  moment  matrix  (given  in  Appendix  2)  is  particularly 
difficult  since  it  involves  five  unrestricted  orbital  summations  plus 
another  two  symmetry  restricted,  orbital  summations  of  two-electron  inte- 
grals. Using  direct  summation  techniques,  the  time  needed  to  construct 
this  matrix  is  roughly  proportional  to  N  .  This  is  a  formidable  computa- 
tional problem,  but  one  that  must  be  accepted  in  favor  of  the  more  manage- 
able matrix  dimensions. 

Fortunately,  the  N  problem  is  not  as  intimidating  as  it  might  seem 
on  first  appearance.  The  "brute-force"  summation  of  two-electron  inte- 
grals in  the  moment  matrices  resums  certain  partial  sums  which  may  appear 
in  more  than  one  term  or  matrix  element.  These  redundant  summations  can 
be  avoided  with  considerable  savings  in  computer  time  by  computing  the 
partial  sums  once  and  reusing  them.  Two  specific  partial  sums  we  have 
employed  are 

[ij|kl]  =     l     <ij|  Iss'xss1  |  |kl>  (2.43) 

s,s' 
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{ij|kl}=     E     <is|  | js'xks*  |  |ls>     .  (2.44) 

S  ,  S  ' 

Since  these  partial  sums  contain  a  double  summation  which  is  performed 
only  once,  the  original  N  problem  is  effectively  reduced  to  N5.  The 
construction  of  the  moment  matrices  is  now  comparable  in  difficulty  to 
the  transformation  of  the  two-electron  integrals  from  the  primitive  basis 

to  the  computational  (usually  Hartree-Fock)  basis  which  is  also  roughly 

5 
proportional  to  N  . 

When  the  number  of  two-electron  integrals  is  too  large  to  be  held 
in  core,  their  random  access  from  peripheral  storage  becomes  relatively 
time  consuming.  The  partial  sums  are  much  more  efficiently  constructed 
from  ordered  lists  of  two-electron  integrals  which  can  be  read  into 
primary  (core)  storage  when  needed.  For  the  partial  sums  defined  above, 
the  two-electron  integrals  must  be  sorted  into  ordered  lists  of  the  type 
<ij|**>  and  <i*|j*>  where  *  indicates  all  orbital  indices  which  yield  a 
non-zero  integral  for  the  corresponding  i,j-th  distribution. 

The  integral  sorts  are  performed  using  the  Yoshimine  sorting  tech- 
nique (Yoshimine,  1973).  Briefly  summarized,  this  technique  involves  a 
partition  of  available  core  into  a  number  of  buffers.  Each  buffer  holds 
integrals  corresponding  to  a  specific  i,j  distribution,  e.g.  <ij|**>. 
(When  the  number  of  distributions  is  large,  several  may  be  held  in  each 
buffer.)  Reading  through  the  two-electron  integral  list,  integrals  are 
then  sorted  into  the  appropriate  buffers.  As  each  buffer  fills,  it  is 
written  to  direct  access,  peripheral  storage  and  assigned  a  record  number. 
All  record  numbers  corresponding  to  integrals  from  the  same  buffer  are 
saved  in  a  "chaining"  array  for  that  buffer.  After  the  entire  integral 
list  has  been  processed  and  all  buffers  have  been  dumped,  it  is  then 
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possible  to  chain  back  through  the  direct  access  records,  copying  inte- 
grals of  the  same  distribution  back  into  core.  These  integrals  may  then 
be  further  sorted  within  distributions,  e.g.  k<l  for  each  i,j,  and 
finally  saved  sequentially  on  a  peripheral  storage  device. 

Diatomic  nitrogen  was  the  first  molecule  to  be  studied  with  the 
[2,1]  Pade'  approximant.  Owing  to  its  abundance  in  the  atmosphere,  nitro- 
gen has  great  chemical  interest  and  has  been  extensively  investigated 
both  experimentally  and  theoretically.   It  is  an  ideal  test  case  for  cal- 
culating ionization  energies  from  a  correlated,  many-electron  formalism 
such  as  propagator  theory  since  both  the  Hartree-Fock  and  AE(SCF)  approx- 
imations incorrectly  predict  the  order  of  the  3a  and  Itt  ionizations 
(Cade  et^  aJL  ,  1966).  Only  when  correlation  corrections  are  included  is 
the  correct  ordering  obtained  (Cederbaum  and  Domcke,  1977  and  references 
therein) . 

A  double  zeta,  contracted  basis  of  Gaussian  type  orbitals  (GTO's) 
was  employed  in  this  calculation.  This  basis  consisted  of  Huzinaga's 
9s, 5p  set  of  primitive  orbitals  (Huzinaga,  1965)  which  was  contracted  to 
4s, 2p  (Dunning,  1970).  This  basis  has  been  optimized  by  Dunning  on  the 
nitrogen  atom  and  is  listed  in  Table  1.  The  corresponding  one-  and  two- 
electron  integrals  were  calculated  at  the  experimental  internuclear 
separation,  R=2.068  a.u.  (Herzberg,  1955),  using  the  MOLECULE  integral 
program  (Almlof,  1974). 

The  Hartree-Fock  calculation  and  the  two-electron  integral  trans- 
formation were  performed  with  the  program  GRNFNC  (Purvis,  1973).  The 
Hartree-Fock  total  energy  with  this  basis  was  E(HF)=  -108.8782  a.u. 
which  is  about  3  eV  higher  than  the  result  of  Cade  et  al.  (1966). 
There  is  also  a  discrepancy  in  the  Hartree-Fock  orbital  energies.  While 
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Table  1.  Contracted  Gaussian  Basis  for  Nitrogen. 

Nitrogen  s  orbital s 

Contraction 
Exponents  Coefficients 

5909.4400  0.002001 

887.4510  0.015310 

204.7490  0.074293 

59.8376  0.253364 

19.9981  0.600576 

2.6860  0.245111 

7.1927  1.000000 

0.7000  1.000000 

0.2133  1.000000 

Nitrogen  p  orbital s 

Contraction 
Exponents  Coefficients 

26.7860  0.018257 

5.9564  0.116407 

1.7074  0.390111 

0.5314  0.637221 

0.1654  1.000000 
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the  calculation  of  Cade  et  aj_.  (incorrectly)  predicted  the  In  orbital 

energy  to  be  0.53  eV  below  the  3a  ,  this  calculation  predicts  the  Itt 

9  u 

energy  to  be  0.05  eV  higher.  The  correct  ordering  of  the  3a  and  Itt 

g     u 
ionizations  with  this  basis  is  merely  fortuitous,  since  based  on  a  total 

energy  criterion,  the  basis  of  Cade  et  al .  is  more  accurate. 

The  next  step  of  the  calculation  involved  the  integral  sorts,  par- 
tial summations,  and  the  construction  of  the  moment  matrices.  The  poles 
and  spectral  density  were  finally  computed  as  outlined  in  the  previous 
section  and  are  presented  along  with  the  [1,0]  results  in  Table  2.  The 
ionization  energies  of  both  approximants  seem  to  be  upper  bounds  to  the 
experimental  results  of  Siegbahn  et  aj_.  (1969),  but  without  exception, 
the  results  of  the  [2,1]  approximant  are  worse  than  the  [1,0]  approximant. 
In  an  attempt  to  incorporate  some  ground  state  correlation  into  the  grand 
canonical  density  operator,  new  occupation  numbers  were  computed  from  the 
spectral  density  and  the  [2,1]  approximant  was  recalculated.  This  cal- 
culation, however,  yielded  no  significant  improvements  in  the  ionization 
energies. 

In  order  to  ascertain  whether  the  poor  results  from  the  [2,1]  approx- 
imant for  nitrogen  are  representative  of  other  calculations  or  just  the 
consequence  of  a  pathological  test  case,  the  water  molecule  was  chosen 
for  a  second  application.  Similarly  to  the  calculation  for  nitrogen,  a 
double  zeta  contracted  basis  of  GTO's  was  also  employed  in  this  calcula- 
tion. Huzinaga's  9s, 5p  primitive  basis  for  oxygen  and  4s  primitive  basis 
for  hydrogen  were  contracted  with  Dunning's  coefficients  to  4s, 2p  and 
2s,  respectively.  The  orbital  exponents  for  the  hydrogen  atoms  were 
scaled  by  1.14  to  more  realistically  represent  the  effective  nuclear 
charge  in  the  molecule, and  the  final  basis  appears  in  Table  3.  Again, 


Table  2.     Principal    Ionization  Energies  for  the  Nitrogen 
Molecule  Resulting  from  the   [1,0]  and   [2,1] 
Propagator  Approximants. 
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Orbital 


la. 


2cj 


3a. 


Itt 


lo 


2a 


[1,0] 

[2,1] 

Exp. 

427.7 

472.4 

409.9 

41.6 

46.5 

37.3 

17.0 

30.4 

15.5 

17.0 

23.1 

16.8 

427.6 

478.8 

409.9 

21.0 

30.9 

18.6 

F)  - 

-108 

8782  H. 

Siegbahn  et  al .  (1969). 
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Table  3.  Contracted  Gaussian  Basis  for  Water. 

Hydrogen  s  sets 

Contraction 
Exponents  Coefficients 

13.3615  0.032828 

2.0133  0.231208 

0.4538  0.817238 

0.1233  1.000000 

Oxygen  s  sets 


Contraction 

Exponents 

Coefficients 

7816.5400 

0.002031 

1175.8200 

0.015436 

273.1880 

0.073771 

81.1696 

0.247606 

27.1836 

0.611832 

3.4136 

0.241205 

9.5322 

1.000000 

0.9398 

1.000000 

0.2846 

1.000000 

Oxygen  p  sets 

Contraction 

Exponents 

Coefficients 

35.1832 

0.019580 

7.9040 

0.124189 

2.3051 

0.394727 

0.7171 

0.627375 

0.2137 

1.000000 
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the  integrals  were  computed  with  the  MOLECULE  program  at  the  equilibrium 
internuclear  geometry,  R(0H)=  1.809  a.u.,  $H0H  =  104.5°  (Benedict  et  al_. , 
1956).  A  total  energy  of  E(HF)=  -76.0082  a.u.  was  computed  with  the 
Hartree-Fock  portion  of  GRNFNC  and  was  followed  by  the  two-electron  inte- 
gral transformation.  Finally,  the  integral  sorts  and  partial  sums  were 
performed,  the  moment  matrices  constructed,  and  the  poles  and  spectral 
density  obtained  for  the  [2,1]  approximant.  The  results  for  both  the 
[1,0]  and  [2,1]  approximants  are  presented  in  Table  4  and  appear  to  be 
upper  bounds  to  the  experimental  ionization  energies.  Once  more,  the 
[2,1]  results  are  consistently  worse  than  the  [1,0]  results.  A  few 
iterations  on  the  occupation  numbers  yielded  no  significant  improvements. 

2.5  Evaluation  of  the  Moment  Conserving  Decoupling 

Formally,  the  moment  conserving  decoupling  is  an  attractive  decou- 
pling procedure.  Being  closely  related  to  the  Pade1  approximant  method, 
this  decoupling  allows  the  application  of  numerous  results  from  the  clas- 
sical moment  problem  to  propagator  theory.   In  particular,  it  was  proven 
that  the  sequence  of  [N,N-1]  approximants  converge  uniformly  to  the 
exact  electron  propagator  in  a  given  basis,  and  it  was  shown  that  these 
approximants  represent  a  variational ly  optimum  choice  of  the  inner  pro- 
jection manifold.  Why  then  are  the  results  of  the  [2,1]  approximant  so 
much  worse  than  the  results  of  the  [1,0]  approximant?  To  answer  this 
question,  it  is  necessary  to  analyze  the  three  approximations  identified 
in  Section  1.3,  namely:  basis  quality,  density  operator,  and  decoupling 
procedure. 

First  of  all,  since  computational  economy  and  not  high  accuracy  was 
the  criterion  for  the  test  calculations  on  nitrogen  and  water,  polarization 
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Table  4.  Principal  Ionization  Energies  for  Water  Resulting 
from  the  [1,0]  and  [2,1]  Propagator  Approximants. 


Orbital 


la. 


[1,01 


2a. 


3a. 


lb, 


lb, 


[2,1] 


E(HF 


Exp.' 


559.4 

619.2 

540.2 

37.0 

44.7 

32.2 

15.4 

29.7 

14.7 

13.8 

32.6 

12.6 

19.5 

29.6 

18.6 

)  =  - 

76 

0082  H. 

Siegbahn  et  al_.  (1969). 
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functions   were    intentionally  excluded  from  the  basis  sets.  Polar- 
ization functions  are  diffuse,  virtual  orbitals  which  can  be  very   impor- 
tant in  describing  electron  relaxation  and  correlation  (Purvis  and  Ohrn, 
1974,  Cederbaum  and  Domcke,  1977).   It  is  reasonable  to  expect  that  the 
addition  of  polarization  functions  will  improve  both  the  [1,0]  and  [2,1] 
approximants  to  varying  degrees;  however,  with  the  same  basis  and  with 
the  same  density  operator,  the  larger  inner  projection  manifold  (if 
judiciously  chosen)  should  yield  a  more  accurate  decoupling.  Since  this 
was  not  the  situation  in  these  test  calculations,  any  improvements  in 
the  basis  sets  did  not  seem  worthwhile. 

Second,  it  is  possible  that  significant  ground  state  correlation 
may  have  been  neglected  with  our  choice  of  the  grand  canonical  density 
operator.  With  the  spin  orbital  annihilation  and  creation  operators 
expanded  in  the  Hartree-Fock  basis  and  using  pure  state  occupation 
numbers  of  zero  or  one,  this  density  operator  yields  the  uncorrected, 
Hartree-Fock  ground  state  average.  Rather  than  explicitly  correlating 
the  density  operator  (as  e.g.  through  perturbation  theory),  an  attempt 
was  made  to  estimate  the  effect  of  correlation  through  the  self-consis- 
tent determination  of  the  occupation  numbers  as  described  in  Section  2.3. 
This  procedure  was  not  pursued  to  true  self-consistency  since  each  iter- 
ation required  a  complete  recalculation  of  the  [2,1]  approximant.   It 
was  obvious,  however,  after  the  first  few  iterations  that  no  significant 
improvements  had  been  obtained. 

Based  on  the  preceding  implications,  the  third  approximation—the 
inner  projection  manifold  truncation—seems  to  be  primarily  responsible 
for  the  poor  numerical  results.  Owing  to  the  complicated  operator  sums 
in  this  manifold,  an  order  analysis  (as  discussed  in  Section  1.6)  is  not 
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readily  possible.  Consequently,  it  is  extremely  difficult  to  identify 
the  problem  with  this  decoupling  procedure.   It  can  only  be  concluded 
that  the  number  of  moments  conserved  is  not  a  useful  criterion  for  decou- 
pling.   This  conclusion  is  consistent  with  the  uniform  convergence  of 
the  [N,N-1]  sequence  since  uniform  convergence  is  not  necessarily  mono- 
tonic,  but  it  suggests  that  more  accurate  decouplings  require  the  incor- 
poration of  more  information  about  the  moment  expansion  than  just  the 
moment  matrices.  The  additional  information  needed  is  indeed  available 
and, in  the  next  chapter,  we  will  demonstrate  how  it  may  be  extracted 
using  perturbation  theory. 


*Babu  and  Ratner  (1972)  reported  the  same  conclusion  which  was  based 
on  an  application  of  their  rational  approximants  to  the  Hubbard  model. 


CHAPTER  3 
DIAGRAM  CONSERVING  DECOUPLINGS 

3. 1  The  Diagrammatic  Expansion  Method 

The  superoperator  formalism  which  is  employed  in  the  previous  two 
chapters  is  by  no  means  the  only  formalism  available  to  formulate  decou- 
pling approximations  for  the  electron  propagator.  Two  commonly  used, 
alternative  methods  are  the  functional  differentiation  method  (see  e.g. 
Csanak  et  al_. ,  1971)  and  the  diagrammatic  expansion  method  (see  e.g. 
Mattuck,  1967,  Fetter  and  Walecka,  1971,  or  Cederbaum  and  Domcke,  1977). 
Of  these  latter  two  methods,  the  diagrammatic  expansion  method  has  proven 
to  be  particularly  effective.  This  method  avoids  some  of  the  algebraic 
tedium  involved  in  deriving  propagator  decoupling  approximations  by 
establishing  certain  rules  for  constructing  and  manipulating  diagrams 
which  represent  the  underlying  algebraic  structure. 

The  diagrammatic  expansion  of  the  electron  propagator  is  usually 
derived  using  time-dependent  perturbation  theory.  The  N-electron  Hamil- 
tonian  is  partitioned  into  an   unperturbed  part  plus  a  time-dependent 
perturbation 

H  =  HQ  +  exp(-e|t|)V  ,  (3.1) 

where  e  is  a  small  positive  quantity.  The  unperturbed  Hamiltonian,  H„, 
is  chosen  to  yield  an  exactly  solvable,  eigenvalue  problem 


w  =  W  '  (3-2: 
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and  the  time  dependence  of  the  unperturbed  eigenstates  is  given  by 

|*0(t)>  =  exp(-iH0t)|<J»0>  .  (3.3) 

In  order  to  simplify  the  remaining  problem  of  finding  the  fully 
perturbed  eigenstates  |¥(t)>,  it  is  convenient  to  introduce  the  "inter- 
action representation"  (Fetter  and  Walecka,  1971)  by  the  transformation 

|Yj(t)>  =  exp(iH0t)|Y(t)>  .  (3.4) 

In  this  representation,  the  Schrbdinger  equation  has  the  form 

i  ft  IVt}>  =  exp(-e|t|)V(t)|Yj(t)>  (3.5) 

where 

V(t)  =  exp(iH0t)Vexp(-iHQt)  .  (3.6) 

The  time  dependence  of  the  interaction  eigenstates  can  be  expressed  as 

|¥j(t)>  =  U£(t,t0)|fI(t0)>  (3.7) 

where  U£(t,t0)  is  the  time-evolution  operator.  Substituting  Eq.  (3.7) 
into  Eq.  (3.5),  the  evolution  operator  is  found  to  satisfy  the  differen- 
tial equation 

i  ft  Ue(t,t0)  =  exp(-e|t|)V(t)Ue(t,t0)  (3.8) 

with  the  initial  condition 

UE(tQ,t0)  -  1  .  (3.9) 

It  is  more  convenient  to  solve  for  U  (t,t„)  by  first  transforming  Eq. 
(3.8)  into  an  integral  equation 
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Ue(t,t0)  =1-1    dtj  exp(-e|t|)V(t1)Ue(t1,t 


O1    ' 


(3.10; 


This  integral  equation  has  the  form  of  the  Volterra  equation  of  the 

second  kind  (Lowdin,  1967)  and  is  solved  iteratively 

t 
Ue(t,t0)  =  1  -  i  J  dtj  exp(-e|t|)V(t1) 


^2 


+  (""•)   j  dtx  J  dt2  exp{-E(|t1|  +  |t2|)}V(t1)V(t2)Ue(t2,t0)      (3. IT 


t0     t0 


00  (  (   1 

1  (-i)n   dt.    dt9 
n=0     >  v  >  l 

l0  l0 


'n-1 


dt  exp{-e( 1 1, |+|tp|+ 


t  I)} 
ni  / 


x  v(t1)v(t2: 


V(t  ) 
v  n7 


tHj^tp^ 


.  >  t 


(3.12) 


Eq.  (3.12)  can  be  generalized  slightly  by  modifying  the  limits  of  inte- 
gration and  introducing  the  time  ordering  operator,  T, 


U  (t,tn)=  Z   (-i)n  -. 
£    °  n=0     n! 


dtl  J  dt2  '  ■  •  J  dtn 
t0     t0 


x  exp{-e(|t1|  +  |t2|+  .  .  .  |tj)}  T[V(t1)V(t2)  .  .  .  V(tn)]  .      (3.13! 


The  time  ordering  operator  rearranges  the  product  of  perturbation  opera- 
tors such  that  the  left-most  term  is  the  latest  in  chronological  order. 

The  perturbed  eigenstates  I ^ j ( tQ)>  can  now  be  expressed  in  terms  of 
the  unperturbed  eigenstates  by  noting  that  as  tQ-*-±°°,  |'Mtn)>^|$n>,  and 
as  tQ  increases  from  -°°  to  zero,  the  perturbation  is  "adiabatically 
switched  on" 
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h<j(o)>  =  ue(o,-°°)|$0>    . 

According  to  a  theorem  of  Gell-Mann  and  Low  (1951),  if 


lim 


KI(0)> 

<*n|U  (o5-«»)|$  >  E  «fn|fT(0) 


Ue(o,-)|*0> 


e+0   *0'  e 


01  'I 


exists,  then  it  is  an  eigenstate  of  H 


Hh'j(O): 


E|Vj(0)> 


'^0I^I(0)>  ~  <$0|fj(0): 


(3.14) 


(3.15) 


:3.16! 


These  results  can  now  be  used  to  determine  the  electron  propagator. 
In  Chapter  1,  the  propagator  was  defined  as  the  ground  state  average  of 
a  time-ordered  product  of  field  operators  in  the  Heisenberg  representa- 
tion 


iG.  .(t) 
ij 


<yH|T[a1(t)at(Q)]|yH> 

<u»  1 5  > 


(3.17; 


Using  Eq.  (3.15)  and  the  fact  that  |TH>= \V, (0)>,  this  average  can  be  ex- 
pressed in  the  interaction  representation  as 


iG.  .(t) 


*n|U J-,t)T[a,(t)at(0)]U  Jt,-»)|*n> 


'01  £ 


r-LL-J  .,--■- 
<$n  U  (-00,00)  $n> 

01  e  '  0 


[3. 18) 


Using  the  expansion  of  the  evolution  operator  (Eq.  (3.13))  and  taking 
the  limit  e-K),  it  can  be  shown  (Fetter  and  Walecka,  1971)  that 


iG..(t)=  M-i)^,  jdtl.  .  .{dtn 


<$0|T[V(t1)  .  .  .  V(tn)a1(t)a](0)]|$0: 


(3.19; 
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The  final  step  in  the  diagrammatic  expansion  method  is  to  expand 
the  numerators  of  each  term  in  Eq.  (3.19)  using  Wick's  theorem  (Wick, 
1950)  and  to  represent  them  diagrammatical ly  (e.g.  Fetter  and  Walecka, 
1971).  The  denominator  of  Eq.  (3.19)  must  also  be  expanded  and  dia- 
grammed, and  when  this  is  done,  all   disconnected    diagrams 
arising  from  the  expansion  of  the  numerator  will  cancel  (Abrikosov 
et  _§_]_. ,  1965). 

Formally,  the  diagrammatic  expansion  method  and  the  superoperator 
formalism  appear  strikingly  dissimilar.  The  diagrammatic  method  is 
formulated  in  the  causal  representation  while  the  superoperator  formal- 
ism utilizes  the  energy  representation.  The  diagrammatic  method  employs 
a  pictorial  representation  of  the  algebraic  structure  while  the  super- 
operator  formalism  emphasizes  the  algebraic  structure  directly.  Yet 
the  primary  goal  of  each  formalism  is  the  same:  an  accurate  prediction 
of  ionization  energies  and  electron  affinities.  Therefore,  the  two 
formalisms  are  inherently  equivalent.   It  is  our  desire  in  this  chapter 
to  explicitly  demonstrate  the  equivalence  between  these  two  formalisms 
and  to  re-examine  the  superoperator  decoupling  approximations  in  terms 
of  a  diagrammatic  analysis. 

3.2  Perturbation  Theory 

The  unifying  feature  of  the  diagrammatic  expansion  method  and  the 
superoperator  formalism  is  perturbation  theory  (Born  and  Ohrn,  1978). 
Since  the  commutator  product  is  distributive  with  respect  to  addition, 
we  can  define  a  partitioning  of  the  superoperator  Hamiltonian  into  an 
unperturbed  part  plus  a  perturbation, 

H  =  HQ  +  V  .  (3.20) 
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One  convenient  partitioning,  which  will  be  shown  to  readily  yield  the 
Hartree-Fock  propagator  as  the  unperturbed  electron  propagator,  is  the 
Mtfller-Plesset  partitioning  (Miller  and  Plesset,  1934).  With  this  par- 
titioning, hL  has  the  form 

H0  =  Z  £rarar  "  ^  E  <rr' |  |  rr'xn  ><n  ,>  (3.21) 

r         r,r' 

and  the  perturbation  is  expressed  as 

v=    E    <rr'}\ss'>[kafal,a    ,a  -5    ,    ,<nr,>a'la]   +  h    I    <rr'\\rr'> 
r,r',s,s'  r  r  s  s  r  s   r   r  s     r^, 

(3.22) 
Of  course,  when  the  commutator  product  is  formed  for  the  superoperators, 
the  constant  term  in  these  definitions  will  cancel. 

Other  parti tionings  of  the  Hamiltonian  may  also  be  assumed  and  may 
lead  to  superior  convergence  properties  (Claverie  et  al . ,  1967).  One 
alternative  partitioning  which  has  been  employed  in  the  perturbation 
calculation  of  correlation  corrections  to  the  total  energy  is  the  Epstein- 
Nesbet  partitioning  (Epstein,  1926,  Mesbet,  1955a,  1955b).  In  propagator 
applications,  the  work  of  Kurtz  and  Ohrn  (1978)  may  be  roughly  interpreted 
in  terms  of  a  partitioning  where  the  unperturbed  Hamiltonian  incorporates 
all  relaxation  contributions  to  the  ionization  energy.  It  is  difficult 
to  define  this  unperturbed  Hamiltonian  explicitly,  but  it  formally  satis- 
fies the  eigenvalue  equation 

H^ak  =  AEk(SCF)ak  (3.23) 

in  contrast  to 

Vk =  ekak  <3-24) 

for  the  M011er-Plesset  partitioning.  The  method  of  Kurtz  and  Ohrn  yields 
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excellent  ionization  energies  and  electron  affinities  with  a  simple 
second-order  self-energy,  however  it  has  not  been  formally  analyzed  in 
detail . 

Corresponding  to  the  partitioning  of  the  superoperator  Hamiltonian, 
we  can  introduce  a  partitioning  of  the  operator  space  defined  by  the 
projection  superoperators  0  and  P, 

6  =  X  |ak)(ak|  =  |a)(a|  (3.25) 

P  =  I  -  0  (3.26) 

These  superoperators  operate  on  elements  of  the  operator  space  through 
the  relations 

OX.  =  E  |ak)(ak|X.)  (3.27) 

PX.  =  X.  -  OX.  (3.28) 

and  are  idempotent  (62  =  6,  P2  =  P),  self-adjoint  (0+  =  6,  Pf  =  P),  and 
mutually  exclusive  (0P=P0=0).  The  superoperator  0  projects  from  an  ar- 
bitrary operator  product  that  part  which  lies  in  the  model  subspace, 
i.e.  that  part  which  is  spanned  by  the  eigenelements  of  f-L.  The  super- 
operator  P  projects  onto  the  orthogonal  complement  of  the  model  subspace, 
i.e.  that  part  which  we  have  no  a  priori  knowledge  about. 

To  obtain  a  perturbation  expansion  of  the  superoperator  resolvent, 
we  consider  its  outer  projection  (Lowdin,  1965)  onto  the  model  subspace, 

G(E)  =  6(EI-H)~16  (3.29) 

=  6(Ei-H0-V)-16  (3.30) 

By  iterating  the  identity 
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(A-B)"1  =   A"1  +  A_1B(A-B)_1  ,  (3.31) 

the  inverse  in  Eq.  (3.30)  can  be  expanded  as 

G(E)  =  (Ei-H0)-10  +  (EI-H0)"10V(EI-H0)"10 

+  (EI-H0)"16v(Ei-H0)"1V(Ei-H0)~16  +  .  .  .  (3.32) 

where  the  property 

[H0,6]_  =  0  (3.33) 

has  been  used.  Now  since  0  plus  P  form  a  resolution  of  the  identity, 
each  resolvent  of  EL  occurring  between  perturbation  superoperators,  V, 
can  be  rewritten  as  a  sum  of  its  projections  on  the  model  subspace  and 
the  orthogonal  complement, 

(EI-Hq)"1  =  (EI-H0)_10  +  (EI-H0)-1P  (3.34) 

=  GQ(E)  +  TQ(E)  .  (3.35) 

With  this  notation,  Eq.  (3.32)  becomes 

G(E)  =  GQ(E)  +  GQ(E)VG0(E)  +  GQ(E)V[GQ(E)  +  TQ(E) ]VGQ(E) 

+  G0(E)V[GQ(E)  +  T0(E)]V[G0(E)  +  TQ(E)]VG0(E)  +  .  .  .       (3.36) 

and  can  be  resummed  to  yield 

G(E)  =  GQ(E)  +  GQ(E)[V  +  VT0(E)V  +  VT0(E)VTQ(E)V  +  .  .  .  ]G(E)    (3.37) 

Defining  the  reduced  resolvent  of  the  full  superoperator  Hamiltonian  as 
T(E)  =  P[a0  +  PCEI-H)?]"-1?  (a^O)  (3.38) 

=  TQ(E)  +  TQ(E)VT(E)  ,  (3.39) 
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Eq.  (3.37)  can  be  written  in  closed  form 

G(E)  *  GQ(E)  +  GQ(E)[V  +  VT(E)V]G(E)  .  (3.40) 

Alternatively,  we  can  define  wave  and  reaction  superoperators  through 
the  equations  (cf.  Lowdin,  1962,  or  Brandow,  1967) 

W(E)  =  I  +  T(E)V  (3.41) 

t(E)  =  VW(E)  .  (3.42) 

The  reduced  resolvent,  wave,  and  reaction  superoperators  introduced 
in  this  section  are  functions  of  the  superoperators  I,  hL,  and  V  and  as 
a  consequence,  operate  in  a  more  complicated  way.  To  apply  a  superoper- 
ator  function  to  an  operator  in  the  operator  space,  it  must  first  be  ex- 
panded in  terms  of  the  superoperators  I,  FL,  and  V  which  are  then  suc- 
cessively applied  to  the  operator.  For  example, 

W(E)X.  =  [I  +  T(E)V]X.  (3.43) 

=  [I  +  Tn(E)V  +  Tn(E)VTn(E)V  +  .  .  .  ]X.  (3.44) 


where 


TQ(E)VXi  =  [E_1I  +  E"2  HQ  +  E"3  H2  +  .  .  .  ]PVX.   ,         (3.45; 


etc. 


3 . 3  Equivalence  of  the  Superoperator  Formalism  and  the  Diagrammatic 
Expansion  Method 

Eqs.  (3.37)  and  (3.40)  represent  the  superoperator  form  of  the  Dyson 

equation  (Dyson,  1949),  and  the  reaction  superoperator  (Eq.  (3.42))  can 

be  identified  as  the  self-energy.  To  demonstrate  that  Eq.  (3.37) 
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corresponds  term  by  term  with  the  diagrammatic  propagator  expansion,  we 
must  first  form  the  operator  average  of  G(E)  to  obtain  the  matrix  Dyson 
equation,  next  evaluate  all  necessary  operator  averages,  and  finally 
diagram  the  resulting  algebraic  formulae.  Owing  to  the  complicated 
operator  averages  that  must  be  evaluated  in  third  and  higher  orders  of 
the  perturbation  superoperator,  the  equivalence  between  these  two  for- 
malisms has  only  been  explicitly  demonstrated  through  third  order  and 
is  assumed  in  all  higher  orders. 

The  matrix  Dyson  equation  is  obtained  by  forming  the  operator  aver- 
age of  G(E)  with  respect  to  the  basis  elements  of  our  model  subspace 

G(E)  =  (a|G(E)a)  (3.46) 

=  G^E)  +  G^E^EJGJE)  ,  (3.47) 

where 

E(E)  =  (aJVa)  +  (a|VTQ(E)Va) 

+  (a|VTQ(E)VT0(E)Va)  +  .  .  .  .  (3.48) 

Since  HQ  was  chosen  to  be  the  Fock  superoperator,  the  appropriate  den- 
sity operator  to  employ  in  the  evaluation  of  the  operator  averages  is 
the  Hartree-Fock  density  operator.  Realizing  that  the  grand  canonical 
density  operator  (Eq.  (1.44))  reduces  to  the  Hartree-Fock  density  oper- 
ator when  pure  state  occupation  numbers  of  zero  or  one  are  chosen,  we 
shall  employ  this  density  operator. 

Beginning  with  the  evaluation  of  matrix  elements  for  the  unperturbed 
propagator,  G^(E),  the  Hartree-Fock  propagator  is  easily  obtained  (cf. 
Section  1.3) 
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G0(E)ij   =    (aj|(EI-H0)-V) 

=  E~l(aJ|a1)  +  E~2(aj  |H0an. )  +  E^a^H^.)  + 
■   f\   +  E-%6^   +  rt?  fi       +   .    .    . 


GQ(E)ir  (E-e^Vj 


(3.49) 
(3.50) 
(3.51) 
(3.52) 


The  evaluation  of  each  term  in  the  self-energy  expansion  requires  the 


initial  evaluation  of  Va . , 


Va.  =  k  Z         <rr'  I I ss ">[a .  ,a  a ' ,a  ,a  ] 

1     r.r'.s.s'        L  i  r  r  s  sJ- 


£   <rs'|  ss'>[a.,a'a  _] 

r,s,s'         ^     r  s 


E      <ir| |ss'>a  a  ,ac  -     I    <is'     ss'xn   ,>a 
r,s,s'  r  s     s       s,s'  s       s 


(3.53) 
(3.54) 


With  this  result,  the  first-order  term  (a.|Va.)  is  obtained  without  much 
additional  effort 
.(1 


E    .  .   =     a.   Va. 

ij  J1      i 


(3.55) 


h    l       <ir| |ss'>Tr{p[a  a c.a_,aT]   } 
r,s,s'  r  s     s     j  + 


•   I     <is'  I  |ss'xn    ,>Tr{p[a    ,a.]    } 
s,s'  J 


(3.56) 


=  7   <ir  js'><n„  ,>5 
r,s  S   rs 


i: 


7   <is ' | I js '><n 


3.57' 


7}l'(E)..    =  0 
ij 


(3.58; 

When  the  effective,  single-particle  potential  used  in  the  unperturbed 
problem  is  the  Hartree-Fock  self-consistent  field  potential,  all  single- 
particle  corrections  vanish  (Bartlett  and  Silver,  1975a). 
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The  evaluation  of  the  second-  and  higher-order  self-energy  matrices 
requires  the  evaluation  of  TQ(E)Vai  and  VTQ(E)Va..  The  first  of  these 
quantities  can  be  expanded  as 

TQ(E)Vai  =  (EI-H0)"1PVai  (3.59) 

=  (EI-H0)"1Vai  -  Z   (EI-H0)"1|ak)(ak|Vai)  (3.60) 

using  Eq.  (3.26).  It  follows  from  the  previous  result  for  (a.|Va.)  that 
the  second  term  in  Eq.  (3.60)  vanishes.  The  first  term  can  now  be  evalu- 
ated by  expanding  the  resolvent  of  hL  and  realizing  that  any  operator 
product  is  an  eigenelement  to  l-L,  i.e. 

Vras'as  =  (er-Es'-es)aras'as  '  <3-61) 

Consequently,  we  obtain 

T0(E)Va.  =  h     T,     ( (E+Er-es-es,)"1<ir|  |ss'>  ajaslas 

-  E  (E-e  )"1<is'  |  |ss'xn  ,>a  (3.62) 

s,s'  b   b 

with  the  help  of  Eq.  (3.45).  The  remaining  application  of  V  and  the 

average  value  evaluation  is  straightforward  and  yields 

£(2)(E)i:J  =  (aj  |VT0(E)Va.)  (3.63) 

_     v  <ir|  I  ss'xss '  I  I  ir>  r,           ,,                                         .              ,_   ... 

=     T,  — t^l — lu —  p-<n  >+i5<n  ><n    ,>-<n   ><n    ,>]              (3.64) 

i    h+c-e-c,  r  s   s     r   s             ' 
r,s,s'   v   r  s  s'  ; 

for  the  matrix  elements  of  the  second-order  self-energy. 

The  Hartree-Fock  average  is  now  obtained  by  choosing  occupation 
numbers  of  zero  and  one.  An  examination  of  the  occupation  number  factor 
in  Eq.  (3.64)  reveals  that  with  this  restriction,  it  will  be  non-vanishing 
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only  when  the  summation  index  r  runs  over  occupied  spin  orbitals  and  s 
and  s'  run  over  unoccupied  spin  orbitals  or  when  r  runs  over  unoccupied 
spin  orbitals  and  s  and  s'  run  over  occupied  spin  orbitals.  Denoting 
a,  b,  c,  .  .  .  as  summation  indices  over  occupied  spin  orbitals;  p,  q, 
r,  .  .  .  for  unoccupied  spin  orbitals;  and  i,  j,  k,  .  .  .  for  unspecified 
spin  orbitals,  Eq.  (3.64)  can  now  be  written  as  two  terms  which  involve 
restricted  spin  orbital  summations 

Z(2)/E)   =  ,,  „   <ia]  |pqxpq|  |ja> 
ij   2a,p,q   ^VVV 

+h     z       <ip||ab><ab||ip>  _  (3_65) 

p,a,b   u  ep~  a~V 

The  conversion  of  Eq.  (3.65)  into  diagrams  is  a  straightforward  pro- 
cedure for  which  we  shall  use  the  rules  and  diagram  convention  of  Brandow 
(1967)  and  Bartlett  and  Silver  (1975b).  This  convention  represents  the 
synthesis  of  the  anti symmetrized  vertices  of  the  Hugenholtz  (1957)  or 
Abrikosov  (1965)  diagrams  with  the  extended  interaction  lines  of  the 
Goldstone  (1957)  diagrams,  and  the  rules  for  constructing  these  diagrams 
are  given  in  Table  5.  The  application  of  these  rules  to  the  terms  in 
Eq.  (3.65)  yields  the  following  diagrams: 

y       <ia||pq><pq||ja>        |  ^  { 
a.p.q   {E+W£q) 

i   v   <ip[  |ab><ab|  |,jp> 

Va,b  t^VVV       <  ^  (3.07; 

These  diagrams  are  precisely  the  same  as  those  obtained  in  the  second- 
order  diagrammatic  expansion  after  a  Fourier  transformation  into  the 
energy  representation  (Cederbaum  and  Domcke,  1977). 
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Table  5.  Rules  for  Constructing  Self-Energy  Diagrams. 

1.  Each  anti symmetrized  two-electron  integral  factor  in  the 
numerator  is  represented  by  an  interaction  line  with  a  vertex 
(dot)  at  both  ends.  The  number  of  interaction  lines  denotes 
the  order  of  the  term. 

2.  Using  the  Dirac  bra-ket  notation,  both  indices  in  the  bra  are 
represented  by  lines  leaving  a  vertex  while  those  of  the  ket 
are  represented  by  lines  entering  a  vertex.  There  must  be  only 
one  outgoing  and  one  incoming  line  per  vertex,  therefore,  assign 
the  index  of  electron  coordinate  one  to  the  left  vertex  and  the 
index  of  electron  coordinate  two  to  the  right  vertex  of  each 
interaction  line. 

3.  Summation  indices  running  over  hole  states  (occupied  orbitals) 
are  directed  downward,  indices  running  over  particle  states 
(unoccupied  orbitals)  are  directed  upward,  and  external  indices 
(not  summed)  are  drawn  horizontally. 

To  Check  Diagrams: 

4.  The  energy  denominator  of  the  diagrammed  expression  should  be 
obtained  by  first  connecting  the  external  lines  and  assigning 
a  factor  of  E  to  this  directed  segment.  Second,  imagine  hori- 
zontal lines  drawn  between  each  pair  of  interaction  lines. 

Each  horizontal  line  corresponds  to  a  multiplicative,  denominator 
factor  obtained  by  summing  the  orbital  energies  of  each  hole 
(downgoing)  line  that  intersects  it  minus  the  sum  of  orbital 
energies  for  particle  (upgoing)  lines  that  intersect  it.  Treat 
the  factor  E  of  the  connected  external  lines  as  an  orbital  energy. 

5.  Numerical  factors  should  be  obtained  by  assigning  a  factor  of  h 
for  each  pair  of  equivalent  internal  lines.  Equivalent  internal 
lines  are  two  lines  which  begin  on  the  same  interaction  line, 
end  on  the  same  interaction  line,  and  go  in  the  same  direction. 

6.  The  overall  sign  factor  should  be  obtained  by  assigning  a  factor 
of  minus  one  to  each  internal  hole  line  segment  and  a  minus  one 
to  each  closed  loop. 
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The  evaluation  of  the  third-order  self-energy  matrix  is  similar  to 
the  second-order  matrix  but  much  more  tedious  and  the  result  is  presented 
in  Appendix  3.  As  was  done  for  the  second-order  expression,  the  occupa- 
tion numbers  must  again  be  restricted  to  zero  and  one  to  obtain  the 
Hartree-Fock  average.  When  this  restriction  is  made,  the  unrestricted 
spin  orbital  summations  in  Appendix  3  will  reduce  to  summations  involving 
occupied,  unoccupied,  and  unspecified  spin  orbitals.  Using  the  algebraic 
identity 


a    _  a    r  i      i  i 

(E-a)(E-b)  "  Ta^bT  [  TI^T  '  T^bT  J 


(3.68; 


it  is  possible  to  combine  terms  in  such  a  way  that  expressions  involving 
only  occupied  and  unoccupied  spin  orbital  summations  are  obtained.  These 
expressions  are  presented  in  Appendix  1.  The  corresponding  diagrams  in 
Appendix  1  again  are  precisely  those  occurring  in  the  third-order,  dia- 
grammatic self-energy  expansion. 


3.4  Diagram  Conserving  Decoupling 

The  wave  and  reaction  superoperators  identified  with  the  help  of 
perturbation  theory  in  Section  3.2  have  special  importance  in  the  develop- 
ment of  decoupling  approximations  for  the  electron  propagator.  As  we 
have  already  seen,  the  reaction  superoperator  generates  the  diagrammatic 
self-energy  expansion.  A  truncation  of  this  expansion  offers  one  viable 
decoupling  scheme.  The  wave  superoperator,  on  the  other  hand,  has  the 
property  of  generating  eigenelements  to  the  full  superoperator  llamilto- 
nian  from  the  eigenelements  of  the  unperturbed  superoperator  Hamiltonian 

(EI-H)W(E)a  =  0  .  (3.69) 
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This  property  is  easily  proven  by  first  using  Eq.  (3.41)  to  expand  W(E) 
and  then  premul tiplying  both  sides  of  Eq.  (3.69)  by  P 

(EI-H)W(E)a  =  (EI-H)a  +  (EI-H)T(E)Va_  (3.70) 

P(EI-H)W(E)a  =  P(EI-H)a  +  P(EI-H)T(E)Va  .  (3.71) 

Using  the  identity 

P(EI-H)T(E)  =  P  (3.72) 

and  the  property  Pa_  =  0,  Eq.  (3.71)  simplifies  to 

P(EI-H)W(E)a  =-  PVa_  +  PVa  =  0  (3.73) 

which  implies  the  validity  of  Eq.  (3.69) 

It  is  of  interest  at  this  point  to  show  a  connection  between  the 
superoperator  formalism  and  the  Equations  of  Motion  (EOM)  method  for 
determining  ionization  energies  (Simons  and  Smith,  1973).   In  this 
method,  one  seeks  solutions  of  the  equation 

[H,Q]_  =  wQ  (3.74) 

which  is  precisely  Eq.  (3.69).  Here  the  operator  Q  is  interpreted  as  a 
correlated  ionization  operator  that  generates,  in  principle,  the  exact 
(N-l)-electron  ion  states  from  the  exact  N-electron  reference  state. 
One  approach  to  solving  Eq.  (3.74)  involves  the  application  of  Rayleigh- 
Schrodinger  perturbation  theory  (Dalgaard  and  Simons,  1977).  By  parti- 
tioning the  Hamiltonian  operator,  expanding  both  the  ionization  operator, 
Q,  and  the  ionization  energy,  to,  in  terms  of  a  perturbation  parameter, 
and  collecting  terms  of  the  same  order,  a  set  of  perturbation  theory 
equations  are  obtained.  The  solution  of  these  equations  yields  an  expan- 
sion for  Q  which  is  analogous  to  the  superoperator  equation 
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h  =  W(E)a  -  (3.75) 

The  only  difference  is  that  E  is  replaced  by  u)Q  which  is  a  consequence 
of  using  Rayleigh-Schrb'dinger  rather  than  Brillioun-Wigner  perturbation 
theory. 

Returning  now  to  the  inner  projection  of  the  superoperator  resolvent, 

G(E)  =  (alhKhKEI-lOhr^hJa)  (3.76) 

we  may  view  Eq.  (3.75)  as  an  alternative  prescription  for  choosing  the 
inner  projection  operator  manifold.  Recalling  from  Section  1.6  that  since 
the  density  operator  describing  the  unperturbed  (model)  problem  does  not 
commute  with  the  full  Hamiltonian,  the  operator  scalar  product  will  not 
in  general  exhibit  Hermitian  symmetry.  Consequently,  we  define 

(h|=(a|Wf(E)  (3.77) 

and  note  that 

(a|Wf(E)  t   (W(E)aJ  .  (3.78) 

Approximate  electron  propagator  decouplings  can  now  be  obtained  by 
truncating  the  expansion  of  the  wave  superoperator, 

W(E)  =  I  +  TQ(E)V  +  TQ(E)VT0(E)V  +  .  .  .  .  (3.79) 

Truncation  of  this  expansion,  with  only  the  superoperator  identity,  triv- 
ially yields  the  Hartree-Fock  propagator,  therefore  we  next  consider 

W(E)  =  I  +  TQ(E)V  .  (3.80) 

Noting  that  the  subspaces  {a.  }  and  (fJfV  =  TQ(E)Va,  }  are  mutually 
orthogonal,  Eq.  (3.76)  can  be  readily  solved  for  G  (E) 
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G_1(E)  =  G^(E)  -  1(E)  ,  (3.81) 

where 

E(E)  =  (aj\/T0(E)Va)(a|VT0(E)(EI-H)T0(E)Va)_1(a[VT0(E)\/a)  (3.82) 
Making  the  following  identifications  from  Section  3.3: 


(a_|VT0(E)Va)  =  EUj(E)  ,  (3.83; 


.(2! 


(a|VTQ(E)(EI-H0)T0(E)Va)  =  Z^;(E)   ,  (3.84) 

and 

(a|VT0(E)VTQ(E)Va)  =  Z(3)(E)   ,  (3.85) 

Eq.  (3.82)  can  be  rewritten 

1(E)  =  E(2)(E)[E(2)(E)  -  E(3,(E)f  V^E)   .  (3.86) 

Expanding  the  inverse  of  Eq.  (3.86),  we  easily  see  that  this  self- 
energy  approximant  coincides  with  the  diagrammatic  expansion  through 
third  order  but  additionally  yields  contributions  to  all  higher  orders. 
If  the  exact  self-energy  is  rewritten  as  a  moment  expansion  in  terms  of 
a  perturbation  parameter,  A, 

X~h(E)   =  E  Ak(a|V(Tn(E)V)ka)  ,  (3.87) 

k=0        u 

we  see  that  Eq.  (3.86)  represents  the  [1,1]  Pade'  approximant  to  this 
expansion.  Owing  to  the  close  connection  between  Pade'  approximants  and 
the  inner  projection  technique  as  demonstrated  in  Chapter  2,  this  result 
is  not  surprising.  These  Pade'  approximants  to  the  self-energy,  however, 
will  have  entirely  different  convergence  properties  than  those  studied  in 
Chapter  2. 
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3.5  Approximations  and  Applications 

Computational  applications  of  the  [1,1]  Pade'  approximant  to  the 
self-energy  require  the  evaluation  of  the  second-  and  third-order  self- 
energy  matrices.  The  second-order  matrix  is  relatively  easy  to  evaluate. 
The  third-order  matrix,  on  the  other  hand,  is  exceedingly  more  difficult 
and  can  presently  be  only  approximately  calculated  without  excessive 
computational  effort.  An  examination  of  the  formulae  in  Appendix  3  re- 
veals that   unlike   the  fourth  moment  matrix  in  the  moment  conserving 
decoupling,  the  third-order  self-energy  matrix  is  energy  dependent. 
This  additional  complication  makes  the  partial  summation  technique  used 
in  the  moment  conserving  decoupling  ineffectual  since  the  third-order 
self-energy  matrix  will  generally  need  to  be  resummed  with  different 
values  of  E  hundreds  of  times  in  the  search  for  poles  of  the  propagator. 

The  first  approximation  that  we  will  examine  is  the  complete  neglect 
of  the  third-order  self-energy  matrix.  With  this  approximation,  the 
[1,1]  approximant  in  Eq.  (3.86)  reduces  to  a  second-order  truncation  of 
the  diagrammatic  self-energy  expansion, 

Z(E)  =  E(2)(E)  .  (3.88) 

This  second-order  self-energy  approximation  is  interesting  not  only  be- 
cause it  contains  the  most  important  relaxation  and  correlation  correc- 
tions to  Koopmans'  theorem  (in  a  perturbation  theoretical  sense),  but 
also  because  it  exhibits  the  same  analytic  form  as  the  exact  self-energy 
(Hedin  and  Lundqvist,  1969,  Cederbaum  and  Domcke,  1977).  Furthermore, 
since  several  second-order,  ionization  energy  calculations  have  been 
reported  in  the  literature,  this  approximation  will  afford  both  a  con- 
venient check  of  new  computer  code  and  the  computational  experience 
necessary  to  implement  more  refined  approximations. 
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The  first  computational  application  of  this  decoupling  approximation 
was  to  the  water  molecule  using  the  same  basis  and  internuclear  geometry 
as  described  in  Section  2.4.  The  results  of  this  calculation  are  pre- 
sented in  Table  6  along  with  the  Koopmans'  theorem,  AE(SCF),  and  experi- 
mental values  for  the  ionization  energies.  Two  ionization  energies  have 
been  tabulated  for  the  2a,  ionization  with  their  corresponding  pole 
strength  (r.  of  Eq.  (1.70))  in  parentheses.  The  occurence  of  two,  rela- 
tively strong  propagator  poles  for  this  ionization  represents  a  breakdown 
in  the  quasi-particle  description  of  inner  valence  ionizations  (Cederbaum, 

1977)  and  makes  assignments  of  principal  and  shake-up  ionizations  ambig- 

* 
uous.   In  general,  the  second-order  ionization  energies  are  quite  en- 
couraging and  represent  significant  improvements  to  each  of  the  Koopmans' 
values.  Furthermore,  these  results  are  comparable  in  accuracy  to  the 
AE(SCF)  results  but  possess  the  convenience  of  being  obtained  in  a  single 
calculation  whereas  the  AE(SCF)  results  required  six  separate  Hartree- 
Fock  calculations. 

The  relatively  poor  agreement  of  the  3a,  and  lb,  ionization  energies 
with  the  experimental  values  in  Table  6  seems  attributable  to  basis  in- 
completeness. Despite  the  lack  of  polarization  functions,  this  suspicion 
is  supported  by  the  facts  that  the  3a,  orbital  is  the  highest  occupied 
orbital  in  that  symmetry  and  that  this  basis  contains  only  two  contracted 
Gaussian  orbitals  of  b,  symmetry.  In  order  to  study  the  basis  dependence 
of  the  second-order  self-energy  approximation,  two  additional  calculations 


*The  ESCA  spectrum  of  the  water  molecule  (Siegbahn  et  aj_.  ,  1969) 
substantiates  this  phenomenon  since  the  2a,  peak  is  quite  broad  and  asym- 
metric. Experimentally,  it  appears  that  trie  lower  energy  ionization 
should  have  a  larger  pole  strength  (in  contrast  with  the  results  of  Table 
5)  since  the  peak  is  skewed  to  higher  binding  energies. 
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Table  6.  Principal  Ionization  Energies  of  Water  Computed 
with  the  14  CGTO  Basis. 


Orbital 

Koopmans 
559.4 

AE(SCF)a 
540.8 

(2) 
£(E) 

Exp.b 

la. 

539.4 

540.2 

z.1 

37.0 

34.6 

34.0  (.61) 
32.6  (.28) 

32.2 

3a, 

15.4 

13.0 

12.9 

14.7 

lbl 

13.8 

11.0 

10.8 

12.6 

lb2 

19.5 

17.8 

18.1 

18.6 

aGoscinski  et  al_.  (1975) 
bSiegbahn  et  al.  (1969). 
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were  performed  with  larger  basis  sets.  The  first  of  these  calculations 
employed  a  26  contracted  orbital  basis  which  augmented  the  original  14 
orbitals  (Table  3)  with  a  set  of  p-orbitals  on  the  hydrogen  atoms  and  a 
set  of  d-orbitals  on  oxygen--all  with  unit  exponents.  The  Hartree-Fock 
total  energy  obtained  with  this  basis  was  E(HF)=  -76.0459  H.  The  second 
calculation  employed  a  38  contracted  orbital  basis  which  included  all  of 
the  orbitals  in  the  26  orbital  basis  plus  an  additional  set  of  diffuse 
p-orbitals  on  the  hydrogen  atoms  (a  =  0.25)  and  a  set  of  diffuse  d-orbit- 
als on  oxygen  (a  =  0.40).  This  basis  yielded  a  Hartree-Fock  total  energy 
of  E(HF)=  -76.0507  H. 

The  most  significant  propagator  poles  calculated  in  the  valence 
region  (0~40  eV)  with  each  of  the  three  water  basis  sets  are  presented  in 
Table  7  along  with  the  second-order  results  of  Cederbaum  (1973a).  The 
inclusion  of  polarization  functions  not  only  improves  the  3a,  and  lb, 
ionization  energies,  it  also  reverses  the  relative  pole  strengths  of  the 
two  dominant  2a^  poles  bringing  the  theoretical  results  into  better  agree- 
ment with  experimental  observations  (see  footnote  on  page  73). 
Cederbaum' s  second-order  results  were  obtained  with  a  basis  comparable 
in  size  and  quality  to  the  26  orbital  basis  in  Table  7.  He  deletes 
several  virtual  orbitals  from  this  basis  before  computing  the  ionization 
energies,  however.  This  approximation  may  account  for  the  small  discrep- 
ancies between  his  results  and  those  reported  here. 

The  formaldehyde  molecule  was  chosen  for  a  second  application  of 
the  second-order  self-energy  approximation.   Ionization  energies  were 
calculated  using  two  basis  sets.  The  first  consisted  of  Huzinaga's  9s, 
5p  primitive  basis  sets  for  oxygen  and  carbon  (Huzinaga,  1965)  contracted 
to  4s  and  2p  functions  with  Dunning's  contraction  coefficients  (Dunning, 
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Table  7.  Basis  Set  Effects  on  the  Ionization  Energies  of  Water 
Computed  with  a  Second-Order  Self-Energy  Approximation, 


Symmetry 

14  CGTO's 

36.5  (.005) 
34.0  (.607) 

32.6  (.279) 
12.9  (.913) 

26  CGTO's 

37.1  (.003) 
33.4  (.288) 
32.1  (.592) 
13.4  (.908) 

38  CGTO's 

33.2  (.231) 
31.9  (.628) 
13.5  (.903) 

Ced.a 

al 

32.9 
13.2 

bl 

34.9  (.005) 
10.8  (.909) 

35.1  (.005) 
11.1  (.904) 

11.2  (.900) 

10.9 

b2 

40.6  (.003) 
18.1  (.931) 

40.8  (.004) 
18.0  (.922) 

18.0  (.919) 

17.7 

E(HF) 

-76.0082 

-76.0459 

-76.0507 

-76.0419 

Cederbaum  (1973a). 
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1970).  The  orbital  exponents  of  Huzinaga's  4s  primitive  basis  for  hydro- 
gen were  scaled  by  a  factor  of  1.2,  and  the  resultant  orbitals  were  con- 
tracted to  2s  functions  as  recommended  by  Dunning.  The  complete  basis 
appears  in  Table  8.  The  second  basis  augmented  the  first  by  the  addition 
of  one  set  of  p-orbitals  on  the  hydrogen  atoms  and  one  set  of  d-orbitals 
on  both  the  oxygen  and  carbon  atoms.  Unit  exponents  were  chosen  for  the 
p-orbitals  on  hydrogen  while  exponents  of  0.8  were  chosen  for  the  d-orbit- 
als. One-  and  two-electron  integrals  were  computed  with  the  MOLECULE 
program  (Almlof,  1974)  at  the  experimental  equilibrium  geometry:  R(C0)= 
2.2825  a.u.,  R(CH)=  2.1090  a.u.,  and  }(HCH)=  116.52°  (Oka,  1960,  Takagi 
and  Oka,  1963),  and  the  Hartree-Fock  calculations  and  two-electron  inte- 
gral transformations  were  performed  with  GRNFNC  (Purvis,  1973).  The 
Hartree-Fock  total  energy  for  the  smaller,  24  orbital  basis  (no  polari- 
zation) was  E(HF)=  -113.8257  H. ,  and  for  the  larger,  42  orbital  basis 
(with  polarization)  E(HF)=  -113.8901  H.  The  Hartree-Fock  orbital  energies 
and  second-order  self-energy  results  for  both  basis  sets  are  presented 
in  Table  9  for  the  principal  ionizations  along  with  the  second-order 
results  of  Cederbaum  et  al_.  (1975)  and  the  experimental  values. 

The  results  in  Table  9  typify  two  general  features  of  ionization 
energy  calculations.  The  first  is  that  Koopmans'  theorem  yields  values 
which  are  usually  higher  than  experimental  ionization  energies.  Second, 
the  inclusion  of  second-order  relaxation  and  correlation  corrections 
generally  overcorrects  the  Koopmans1  estimate  and  yields  values  which 
are  usually  lower  than  experiment.  For  several  ionizations  in  Table  9, 
the  second-order  deviations  from  experiment  are  as  large  as  the  Koopmans' 
values  only  opposite  in  sign.-  Although  it  is  possible  that  the  larger, 
polarized  basis  used  in  the  second  calculation  may  still  lack  adequate 
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Table  8. 

Contracted  Gaussian  Basis  for  Formaldehyde. 

Carbon 

s  sets 

Oxygen 

s  sets 

Contraction 

Contraction 

Exponents 

Coefficients 

Exponents 

Coefficients 

4232.6100 

0.002029 

7816.5400 

0.002031 

634.8820 

0.015535 

1175.8200 

0.015436 

146.0970 

0.075411 

273.1880 

0.073771 

42.4974 

0.257121 

81.1696 

0.247606 

14.1892 

0.596555 

27.1836 

0.611832 

1.9666 

0.242517 

3.4136 

0.241205 

5.1477 

1.000000 

9.5322 

1.000000 

0.4962 

1.000000 

0.9398 

1.000000 

0.1533 

1.000000 

0.2846 

1.000000 

Carbon 

p  sets 

Oxygen 

p  sets 

Contraction 

Contraction 

Exponents 

Coefficients 

Exponents 

Coefficients 

18.1557 

0.018534 

35.1832 

0.019580 

3.9864 

0.115442 

7.9040 

0.124189 

1.1429 

0.386206 

2.3051 

0.394727 

0.3594 

0.640089 

0.7171 

0.627375 

0.1146 

1.000000 

0.2137 

1.000000 

Hydrogen  s  sets 


Exponents 

Contraction 
Coefficients 

19.2406 
2.8992 
0.6535 
0.1776 

0.032828 
0.231208 
0.817238 
1.000000 
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Table  9.  Principal  Ionization  Energies  for  Formaldehyde. 
24  42 


Orbital 

Koopmans 
560.12 

(2) 
£(E) 

Koopmans 

(2) 
£(E) 

Ced.a 

Exp. 

la  j 

538.93 

559.81 

538.62 

- 

539. 43b 

2a, 

309.09 

297.27 

308.87 

296.90 

- 

294. 21b 

3a, 

38.94 

33.66 

38.18 

32.56 

- 

34. 2C 

4a  j 

23.39 

20.97 

23.30 

21.03 

- 

21.15C 

5a, 

17.29 

13.98 

17.38 

14.38 

14.42 

16. 2d 

lbl 

14.56 

13.83 

14.45 

13.72 

13.50 

14. 5d 

lb2 

19.47 

17.16 

19.08 

17.07 

16.63 

17. 0d 

2b2 

12.06 

9.04 

11.93 

9.30 

9.25 

10. 9d 

E(HF) 

-113 

8257 

-113 

8901 

-113.9012 

Second  order  results  of  Cederbaum  et  al_.  (1975). 

bJolly  and  Schaaf  (1976). 

cHood  et  al_.  (1976). 

Estimated  center  of  gravity  (Cederbaum  and  Domcke,  1977)  from  spectrum 
of  Turner  et  al.  (1970)  . 


polarization  functions,  the  rather  large  discrepancies  between  the 
second-order  results  and  experiment  more  probably  indicate  that  third- 
(and  higher)  order  self-energy  corrections  are  now  sizable.  The  general 
conclusion  that  a  second-order  self-energy  approximation  is  inadequate 
for  an  accurate  calculation  of  ionization  energies  has  been  previously 
concluded  by  Cederbaum  (1973b)  and  necessitates  a  re-examination  of  the 
approximation  made  in  Eq.  (3.88). 

Rather  than  completely  neglecting  the  third-order  self-energy  matrix, 
let  us  now  consider  an  approximation  that  includes  at  least  part  of  these 
contributions.  Which  third-order  self-energy  diagrams  should  be  included? 
There  are  two  well-established  results  that  are   relevant  to  this  ques- 
tion: Studies  of  the  electron  gas  model  have  shown  that  in  the  limit  of 
high  electron  density,  the  so-called  ring  diagrams  dominate  the  self- 
energy  expansion  (Pines,  1961),  while  in  the  limit  of  low  electron  den- 
sity, the  so-called  ladder  diagrams  dominate  (Galitskii,  1958).   In  order 
to  determine  whether  atomic  and  molecular  self-energies  can  be  approxi- 
mated by  specific  third-order  diagrams  (e.g.  rings  or  ladders),  we  need 
to  evaluate  all  third-order  diagrams  for  some  representative  systems. 
Cederbaum  (1975)  has  done  this  for  several  simple  systems  and  has  found 
that  both  ring  ajid  ladder  diagrams  dominate  the  third-order  self-energy. 
This  result  implies  that  atoms  and  molecules  lie  somewhere  between  the 
high  and  low  density  extremes.   It  is  therefore  essential  to  include 
both  ring  and  ladder  diagrams  in  any  third-order  self-energy  approxima- 
tion. These  diagrams  are 

rings  ladders 

"L  71 


..3  1-3 


3.89 
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and  correspond  to  the  algebraic  expressions  labeled  A,  B,  C,  and  D  in 
Appendix  1. 

We  include  six  additional  diagrams  in  our  third-order  self-energy 
approximation  because  of  the  computational  efficiency  with  which  they 
are  evaluated.  These  diagrams  are  the  energy  independent  diagrams 


*ml:ih   ,,,, 


corresponding  to  expressions  M-R  in  Appendix  1.  For  these  six  diagrams, 
it  is  feasible  to  employ  the  partial  summation  technique  since  they  must 
be  evaluated  only  once. 

Approximating  the  full  third-order  self-energy  matrix  by  only  ring, 
ladder,  and  constant  energy  diagrams,  let  us  now  consider  the  solution 
of  the  Dyson  equation  with  the  [1,1]  Pade'  approximant  to  the  self-energy 
expansion.  Owing  to  the  fact  that  the  inner  projection  manifold  from 
which  the  [1,1]  approximant  was  derived  is  energy  dependent  (Eq.  (3.80)), 
the  simple  analytic  form  of  the  self-energy  eigenvalues,  illustrated  in 
Fig.  2,  is  lost.  Furthermore,  the  self-energy  poles  are  now  given  by 

det  (E(2)(E)  -  E(3)(E))  =  0  (3.91) 

rather  than  by  an  eigenvalue  problem  and  are  consequently  more  difficult 
to  obtain.  For  these  reasons,  the  pole  search  described  in  Chapter  1 
and  used  with  the  second-order  self-energy  approximation  is  no  longer 
an  efficient  or  reliable  procedure.  An  alternative  method  of  solution 
used  in  the  following  applications  was  to  use  the  Hartree-Fock  orbital 
energy  as  an  initial  guess  to  the  propagator  pole  and  to  iterate  Eq. 
(1.67)  to  convergence.  When  convergent,  this  procedure  invariably  yields 
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a  principal  propagator  pole  and  its  corresponding  pole  strength.  Al- 
though the  [1,1]  self-energy  approximant  does  not  quarantee  a  positive 
pole  strength,  this  was  never  a  problem  in  any  of  the  calculations  re- 
ported here. 

The  principal  ionization  energies  for  the  water  molecule  were  cal- 
culated using  both  the  14  and  26  CGTO  basis  sets  in  order  to  evaluate 
the  [1,1]  Pade'  approximant  to  the  self-energy  expansion,  and  the  results 
appear  in  Table  10.  The  most  significant  feature  of  these  results  is  that 
each  ionization  energy  has  been  shifted  from  its  second  order  value  to 
higher  energy  and  is  now  in  better  agreement  with  the  experimental  value. 
It  is  further  noticed  that  the  valence  ionization  energies  are  still 
smaller  than  the  experimental  values  while  the  la1  (core)  ionization 
energy  is  now  larger  than  experiment.  Apparently,  the  diagrams  included 
in  the  third-order  self-energy  matrix  overestimate  the  actual  relaxation 
and  correlation  effects  for  this  ionization. 

Convergence  difficulty  was  experienced  for  the  2a.,  ionization 
energy  using  the  14  orbital  basis.  A  schematic  plot  of  VL      (E)  is  pre- 
sented  in  Fig.  3  and  reveals  that  there  are  no  propagator  poles  in  this 
energy  region.  This  anomaly  is  no  doubt  a  consequence  of  some  quirk  in 
the  basis  since  the  26  orbital  basis  yields  a  very   accurate  2a,  ioniza- 
tion energy. 

3 . 6  Evaluation  of  the  Diagram  Conserving  Decoupling 

The  algebraic  structure  of  the  superoperator  formalism  has  been 
successfully  exploited  in  this  chapter  to  yield  several  new  insights  into 
the  decoupling  problem.  The  application  of  perturbation  theory  has 
demonstrated  that  the  electron  propagator  equation  of  motion  can  be 
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Table  10.  Comparison  of  Principal  Ionization  Energies  for  Water 

Obtained  with  the  Second-Order  and  the  [1,1]  Self-Energies 
Using  the  14  and  26  CGTO  Basis  Sets. 


14 


26 


Orbital 


la. 


2a. 


3a. 


lb, 


lb. 


(2) 
E(E) 

[1,1] 

(2) 
Z(E) 

[1,1] 

Exp.a 

539.4 

541.6 

539.2 

540.9 

540.2 

34.0  (.607) 
32.6  (.279) 

no  convergence 
(see  text) 

32.1 

32.2 

32.2 

12.9 

13.4 

13.4 

13.6 

14.7 

10.8 

11.1 

11.1 

11.3 

12.6 

18.1 

18.4 

18.0 

no  results 

18.6 

Siegbahn  et  al.  (1969) 


Figure  3.  A  sketch  of  W2a,  in  the  energy  region  of  the 
2a j  ionization  obtained  with  the  [1,1]  self- 
energy  approximant  using  the  14  CGTO  basis. 
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W2a(E) 


-1.215         -1.175 


resummed  to  yield  the  equivalent  of  the  diagrammatic  expansion.  This 
resummation  also  allows  the  identification  of  wave  and  reaction  super- 
operators  which  have  special  importance  in  decoupling  approximations. 
We  have  shown  that  when  the  inner  projection  manifold  of  the  superoper- 
ator  resolvent  is  chosen  to  consist  of  the  first-order  truncation  of  the 
wave  superoperator,  a  [1,1]  Pade'  approximant  to  the  self-energy  expan- 
sion is  obtained.  This  approximant  is  correct  through  third  order  and 
contains  a  geometric  approximation  to  all  higher  orders.   In  general, 
the  Nth-order  truncation  of  the  wave  superoperator  will  yield  an  [N,N] 
Pade'  approximant  which  is  correct  through  the  (2N+l)st  order  in  the 
self-energy  expansion.  One  final  insight  afforded  by  this  decoupling 
is  the  realization  that  electron  correlation  can  be  described  exclusively 
in  the  operator  space.  We  argued  in  Chapter  1  that  when  the  propagator 
was  defined  as  a  single-time  Green's  function,  the  density  operator  was 
arbitrary.  We  have  now  demonstrated  in  this  chapter  that  any  desired 
order  in  the  self-energy  expansion  may  be  obtained  using  as  a  specific 
choice, the  uncorrelated,  Hartree-Fock  density  operator. 

Computational  applications  of  the  diagram  conserving  decoupling 
have  been  encouraging.  These  applications  have  confirmed  previous  con- 
clusions (Cederbaum,  1973b)  that  a  second-order  self-energy  is  generally 
inadequate  for  obtaining  accurate  ionization  energies.  It  is  important 
if  not  essential  that  third-order  ring  and  ladder  diagrams  be  included 
in  any  self-energy  approximation  although  the  errors  arising  from  basis 
incompleteness  may  be  of  equal  magnitude  and  hence  cannot  be  ignored. 
The  inclusion  of  the  third-order  ring,  ladder,  and  constant  energy  dia- 
grams in  the  [1,1]  self-energy  approximant  has  succeeded  in  improving  the 
second-order  results  but  even  these  results  are  not  consistently  better 
than  the  Koopmans'  theorem  values. 
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One  important  feature  of  the  [1,1]  self-energy  approximant  is  that 
even  though  it  is  constructed  from  only  the  second-  and  third-order  self- 
energy  matrices,  it  contains  a  geometric  approximation  of  all  higher 
orders  in  the  self-energy  expansion.  Certainly,  some  fourth-  or  higher- 
order  terms  may  be  just  as  important  as  third-order  terms;  therefore, 
this  approximation  is  highly  desirable.  The  fourth-  and  higher-order 
terms  arising  from  the  [1,1]  self-energy  approximant,  however,  are   not 
readily  analyzed  diagrammatical ly.  In  fact,  being  a  purely  algebraic 
approximation,  the  [1,1]  approximant  may  not  yield  any  valid  fourth-  or 
higher-order  diagrams.  Given  the  fact  that  ring  and  ladder  diagrams 
dominate  the  third-order  self-energy  matrix,  one  can  argue  that  they  may 
also  dominate  the  higher  orders  of  the  self-energy  expansion.  An  appro- 
priate modification  of  this  decoupling  scheme  might  then  allow  the  sum- 
mation of  these  specific  diagrams  in  all  orders.  Approximations  of  this 
type  are  referred  to  as  renormalized  decouplings  and  are  examined  within 
the  superoperator  formalism  in  the  next  chapter. 


CHAPTER  4 
RENORMALIZED  DECOUPLINGS 

4. 1  Renormal ization  Theory 

In  Chapter  3,  we  tacitly  assumed  that  the  application  of  perturba- 
tion theory  to  the  calculation  of  ionization  energies  and  electron 
affinities  was  valid  and  that  the  resulting  self-energy  expansion  was 
convergent.  Historically  however,  it  was  discovered  that  in  both  the 
nuclear  many-body  problem  and  the  electron  gas  model,  the  simple  self- 
energy  expansions  are  divergent.  In  order  to  remove  these  divergencies, 
it  is  necessary  to  sum  certain  appropriate  classes  of  diagrams  to  all 
orders.  This  method  of  partial  summations  is  known  as  renormal ization 
theory  (see  e.g.  Kumar,  1962  or  Mattuck,  1967)  and  may  be  viewed  as  an 
analytic  continuation  of  the  perturbation  expansion.  Although  a  variety 
of  renormal ization  procedures  exist,  such  as  propagator  renormal izations, 
interaction  renormal izations,  and  vertex  renormal izations ,  the  distinc- 
tions mainly  depend  on  the  types  of  diagrams  included  in  the  partial 
summation  and  are  not  particularly  important  for  our  consideration. 

One  renormal ization  that  we  are  already  familiar  with  is  the  [1,1] 
self-energy  approximant  derived  in  the  preceding  chapter.   In  fact,  any 
rational  self-energy  approximant  may  be  regarded  as  a  renormal ization 
since  its  geometric  expansion  will  approximate  all  orders  of  the  pertur- 
bation expansion.  One  problem  encountered  with  the  [1,1]  approximant 
and  that  occurs  in  general  for  rational  approximants  derived  via  purely 
algebraic  considerations  is  that  their  geometric  expansions  may  contain 


no  readily  identifiable  diagrams  (at  least  beyond  the  lowest  orders). 
Since  specific  diagrams  often  dominate  the  self-energy  expansion  (such 
as  ring  and  ladder  diagrams  for  atoms  and  molecules)  it  is  valuable  to 
investigate  whether  the  superoperator  formalism  can  be  adapted  to  yield 
renormalized  self-energy  expressions  that  sum  specific  diagrams.  The 
solution  as  we  shall  see  is  rather  simple  and  involves  a  restriction  in 
the  types  of  operator  products  allowed  to  span  the  orthogonal  complement 
of  the  model  subspace.  As  a  specific  example,  the  two  particle-one  hole 
Tamm-Dancoff  approximation  (2p-h  TDA),  (Schuck  et  al_. ,  1973,  Schirmer 
and  Cederbaum,  1978),  is  derived  from  an  effective  interaction  which  is 
logically  obtained  by  a  projection  of  the  perturbation  superoperator  onto 
the  subspace  spanned  by  2p-h  type  operators  (Born  and  Ohrn,  1979). 
Finally,  the  diagonal  approximation  to  the  full  2p-h  TDA  self-energy 
previously  derived  and  applied  to  the  calculation  of  ionization  energies 
is  shown  to  neglect  terms  which,  in  fact,  are   diagonal  and  are  necessary 
to  prevent  an  overcounting  of  all  diagrams  containing  diagonal  ladder 
parts. 

4.2  Derivation  of  the  2p-h  TDA  and  Diagonal  2p-h  TDA  Equations 

Recalling  some  of  the  results  of  the  previous  chapter,  we  had  ob- 
tained the  matrix  Dyson  equation 

G(E)  =  G^E)  +  Gq(E)Z(E)G(E)  ,  (4.1) 

where  the  self-energy  matrix,  EJE),  had  the  following  expansion 

Z(E)  =  (a|V+VT0(E)V+VT0(E)VT0(E)V+  .  .  .  |a)  .  (4.2) 
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Introducing  the  reduced  resolvent  of  the  full  superoperator  Hamil tonian, 
T(E),  which  is  just  a  projection  of  the  superoperator  resolvent  on  the 
orthogonal  complement 

T(E)  =  P[a6+(EI-H0)P-PVP]"1P  ,  (4.3) 

the  self-energy  expansion  was  written  in  closed  form 

2(E)  =  (a|Va)  +  (aJvT(E)Va)  .  (4.4) 

It  was  further  shown  that  when  the  grand  canonical  density  operator  is 
used  to  evaluate  the  operator  averages,  the  first-order  term  vanishes. 
When  P  is  the  exact  projector  of  the  orthogonal  complement, 

P  =  I  -  0  ,  (4.5) 

the  term  PVP  in  Eq.  (4.3)  is  responsible  for  generating  the  operator 
products  that  span  this  subspace.  The  expansion  of  this  term  from  the 
inverse  and  its  repeated  application  in  each  order  of  the  perturbation 
expansion  yields  larger  and  larger  operator  products  which  are  only 
limited  by  the  number  of  electrons  in  the  reference  state.   If  instead 
of  allowing  all  possible  operator  products,  we  restrict  them  to  some 
simple  types  which  occur  in  each  order,  it  may  be  possible  to  identify 
and  sum  specific  diagrams  in  all  orders  of  the  perturbation  expansion. 

The  restriction  of  the  operator  products  in  the  orthogonal  comple- 
ment is  achieved  by  approximating  the  orthogonal  projector  as 

P  =  ll)(ll  (4.6) 

where  the  manifold  (f)  contains  the  desired  operator  products.  The  pro- 
jector P  now  has  the  effect  of  projecting  from  the  perturbation  expansion 
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only  those  operator  products  which  lie  in  the  subspace  spanned  by  {f}. 
The  approximation  to  P  in  Eq.  (4.6)  must  of  course  preserve  the  proper- 
ties of  the  exact  projector  and  should  be  idempotent,  self-adjoint,  and 
orthogonal  to  0. 

Our  previous  experience  with  the  operator  product  decouplings 

suggests  that  the  set  of  triple  operator  products  {a.^a  }  be  chosen 

k  1  m 

as  a  first  approximation  to  P.  There  is  a  stronger  motivation  for  using 
this  operator  product,  however.  If  the  third-order  ring  and  ladder  dia- 
grams in  Eq.  (3.70)  are  examined,  it  can  be  seen  that  between  any  two 
interaction  lines  there  occurs  only  two  particle  lines  (upgoing)  and  one 
hole  line  (downgoing)  or  vice  versa.  This  implies  that  the  intermediate 
or  virtual  states  that  are  represented  by  these  diagrams  consist  of  only 
2p-h  or  2h-p  excitations  of  the  reference  state.  Both  of  these  excita- 
tions are  described  with  the  triple-operator  products. 

The  set  of  triple  products  {a.a,a  }  is  not  orthogonal  to  the  simple 
operators  of  the  model  subspace,  hence  these  two  subspaces  must  be  or- 
thogonal ized.  Using  the  Gram-Schmidt  orthogonal ization  procedure  (see 
e.g.  Pilar,  1968),  we  define 

_!.     .1.  J. 

fnm  =  Ni,i2  [a,a,a  -S(a    la,  a, a   )a   ]  (4.7) 

klm         klnr   k   1   m         n1    k   1   nr    n  v        ' 

=  H.J   [a,  an a  +6.    <n.  >a,-6.  -,<n,  >a   1  (4.8) 

klm1   k  1   m     km     k     1     kl     k     m  v       ' 

where 

NUm  =  <n.  >-<n.  ><n1>-<n,  ><n  >+<n,xn  >   .  (4.9) 

klm  k         k       1  k       m         1       m  v        ' 

The  projector 

P=kj5Jfklm)(fklJ        ,      Km  (4.10) 
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is  now  idempotent,  self-adjoint,  and  orthogonal  to  0.  The  projection  of 
the  perturbation  superoperator  on  the  2p-h  subspace,  PVP,  which  occurs  in 
Eq.  (4.3)  can  now  be  regarded  as  an  effective  interaction.  The  expansion 
of  Eq.  (4.3)  with  P  defined  as  in  Eq.  (4.10)  should  yield  all  diagrams 
containing  2p-h  and  2h-p  excitations  of  the  reference  state. 

The  necessary  operator  averages  needed  to  evaluate  PVP  are: 

(a'1a1Iaml|5a])  =  (a,  IVa+.a^.)*  =  N^  v<lk'  |  |«n' V>     (4.11) 

and 

(ak'al  *am'  'Vakalam)  =  Nk'l'm'{<mlMm,1,>«kk.(l-<V<nl>) 

-<k'ni|  I  km '  >6,  -,  ,  (<n,  >-<n  >)-<k'l-|  Ikl  ">6  ,(<n,>-<n,>) 
11     11    k    m      ' '     mm  v  k    1  ; 

+<k 'ml  I  kl  ' >6  ,1(<n.>-<n  >)+<k'l  I  I  km'>6  -, ,  (<n,  >-<n,>) 

11    m  1   k    m      '  '     ml    k    1  ' 

+<k'm|  |1  'm'>6kl<n  >+<k'l  |  |m'T>6,  <n,  >} 

+Nklm{<ml  |  |kT>6ni,k,<nni,>+<ml  |  |m'k>6]  .k.<n]  ,>}  (4.12) 

Substituting  these  expressions  and  performing  some  cancelation  yields: 

p^\jimk,j,jm/(klNk'rm'<™'iln.'r>6kk,(i-<„ni>-<ni>)|fk,,,m,)(fk1m| 

1  <  m  1  '<  m 

^j,mk^^,m'N»mn^'m,<k'm|ikm'>6'1'(<v"<v),fk•',m,)(fk,J 

1  <  m  1 '  <  m ' 

-     £  >'■         N^N?.-,.    .<k,l||kr>6      ,(<n.  >-<n,>)|f,  ,,  ,    ,)(f.,    I 

k   1   m  k'    1'   m'  k'l'm'  M  mm'v      k         1    "    k   1   m         klm1 

1   <  m     1  '  <  m ' 

+     T.  I         N.-'-fN;2,,,    ,<k'm|  |kl'>6n    ,  (<n.  >-<n  >)  If.  ,, ,    , )  (f. ,    I 

klmk'l'm'  m 

1   <  m     1  '  <  m ' 

+     z  l        NMmNb-Tm-<k,1M'<rnl>6  11(<n.>-<n1>)|f,  ,,,    ,)(f.,    I 

klmk'l'm' 
1    <  m     l'<  m'  (4.13) 
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Additional  simplification  can  be  achieved  at  this  point  by  anti -commuting 
the  operators  a] ,  and  am,  in  I ^^ .  -j . m . )  of  the  last  two  terms  in  Eq.  (4.13! 
with  the  appropriate  change  of  sign, 

IVl'm')  =  "  lfkVl'}  •  <4-14) 

After  interchanging  dummy  indices  1  '<->-  m',  the  fourth  and  fifth  terms  be- 
come equal  to  the  second  and  third  terms,  respectively,  and  since  the 
diagonal  terms  1'=  m'  vanish,  the  summations  with  l'<  m1  and  l'>  m'  can 
be  combined  as  unrestricted  summations  over  1'  and  m'.  Since  the  first 
term  in  Eq.  (4.13)  is  obviously  symmetric  in  1'  and  m',  the  restriction 
l'<  m'  in  that  term  may  be  removed  by  multiplying  the  sum  by  a  factor 
of  k-     The  remaining  restriction,  1  <m,  may  also  be  removed  by  introduc- 
ing another  factor  of  k   since 


1-.   V, 


PVP  =  1  I  Nkl2mNk1l'm,{1'2<mlllm'1'>6kk'(1-<nm>-<nl>) 

k,l,m  k',r,m'  Klm  K  '  m  kk     m    ' 

1  <  m 

-<k'm| |km'>611 ,(<n,  >-<nm>)-<k'l I Ikl '>6  ,(<n,  >-<n1>)}|f,  M ,  ,)(f.,  I 
11    k    m      ' '    mm  v  k    1  '  '  k' lm  v  klm1 

(4.15) 

is  symmetric  in  these  indices  as  can  be  verified  by  interchanging  Wm 

and  relabeling  dummy  indices  l'-w-m'.  Expanding  the  ket-bra  superoperator 

hi  E    lfk'Tm')(fklml  (4-16^ 

k,l,m  kM',m'  K  '  m    klm 

out  of  the  inverse,  evaluating  the  remaining  operator  averages,  and 

resumming  the  expansion  yields  the  2p-h  TDA  self-energy 

2p-h  TDA  ,   , 

2(E).  ■   =  h     Y,  l         N,s,  N,3,, ,  ,<ik|  |lm> 

1J     k,l,m  k',l',m'  klm  kVm 

x  {(fKEi-H^fJ-dlVDJ-^^^^^Tm'llJk^  (4.17) 

where 
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(fk'Tm'KEI-Ho)W=(E+ek-el-eJ6kk'6ir6mm-  ■ 

f^H'm'l5WBNi;W'rm'^ralllm,1,>fikk'<1-<nm>-<nl>) 

-<k'm|  |km,>611,(<nk>-<n|Ti>)-<k*1|  |kl  ' >6mm ,  (<n|<>-<n-| > ) }      (4.18) 

Although  different  in  appearance,  this  self-energy  expression  is  formally 
the  same  as  that  obtained  by  Purvis  and  Ohrn  (1975a)  using  the  operator 
product  decoupling  and  by  Cederbaum  (1975)  and  Schirmer  and  Cederbaum 
(1978)  using  the  diagrammatic  method.  The  present  derivation  clearly 
illuminates  the  parallelism  between  the  two  formalisms. 

Owing  to  the  large  dimension  of  the  {f,  ,  }  operator  subspace  and 
the  associated  difficulty  in  diagonal izing  (f|Vf),  computational  appli- 
cations of  the  2p-h  TDA  have  usually  involved  additional  approximations. 
One  approximation  which  has  facilitated  computational  applications  is 
known  alternatively  as  the  shifted  Born  collision  (SBC)  approximation 
(Purvis  and  Ohrn,  1974,  1975a)  or  the  diagonal  2p-h  TDA  (Cederbaum, 
1974,  1975  and  Cederbaum  and  Domcke,  1977).  This  purportedly  "diagonal" 
approximation  restricts  the  spin-orbital  summation  indices  in  Eq .  (4.13) 
to  k'  =  k,  l'  =  l,  and  m'=m  thereby  neglecting  the  last  two  summations 
and  yielding  the  following  self-energy  expression 

2p-h  TDA  ..  |  ,,   ,  ,,., 

E(E)..  -H     I       N.,  <;k||1m><H  jk>  (4   } 

1J        k,i,m  k1m  (£+v£rv-A 

where 

A=<ml | |ml >( l-<n  >-<n, >)-<km| |km>(<n, >-<n  >)-<kl | | kl > (<n. >-<n, >) .      (4.20) 
in  i  K  ill  K  1 

By  neglecting  the  last  two  summations  in  Eq.  (4.13),  however,  this  approx- 
imation actually  neglects  some  diagonal  contributions  to  (f.  ,-,,  ,  |Vf, ,  ). 
3  3  v  k  1  m  '  klm; 

As  we  have  explicitly  demonstrated  in  the  derivation  of  Eq.  (4.17)  and 
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(4.18),  the  2p-h  TDA  self-energy  sums  are  symmetric  in  both  1 ,m  and 
l',m';  consequently,  the  last  two  summations  in  Eq.  (4.13)  contain  pre- 
cisely the  same  contributions  as  the  second  and  third  summations,  re- 
spectively. If  the  diagonal  approximation  (k'=k,  l'=l,  m'=m)  is  made  in 
Eqs.  (4.17)  and  (4.18),  this  symmetry  is  properly  accounted  for  and  the 
resulting  self-energy  expression  is 

2p-h  TDA  .,  ,  ,,   ,  I,, 

Z(E)     =  h     T     N    <ikl[1m><1m  |jk>  (4  n) 

{)"        \j>^I^WJ^  (4'21) 

where 

A=?2<ml  I  I ml >( l-<n  >-<n,>)-<km| I km>(<n, >-<n  >)-<kl I |kl>(<n, >-<n,>) .  (4.22) 

II    \    m    i  )  II    v  y  m  /     II    v  y  ]  /   \     / 

Eq.  (4.22)  differs  from  Eq.  (4.20)  by  the  factor  of  h   in  the  first  term. 
The  inclusion  of  this  factor  in  the  diagonal  2p-h  TDA  is  necessary  to 
prevent  an  overcounting  of  diagrams  with  diagonal  ladder  parts  (see  next 
section)  and  typically  shifts  ionization  energies  0.3-0.4  eV  higher  in 
energy  (Born  and  Ohrn,  1979). 

4.3  Diagrammatic  Analysis 

In  order  to  determine  precisely  which  self-energy  diagrams  are   in- 
cluded in  the  2p-h  TDA  self-energy,  it  is  necessary  to  expand  the  effec- 
tive interaction  matrix,  (jfJVfj,  from  the  inverse  in  Eq.  (4.17)  and  dia- 
gram the  resulting  algebraic  expressions  in  each  order.  The  expansion 
of  Eq.  (4.17)  yields  the  following  terms  in  lowest  orders 

2p-h  TDA  , 

E(E).  ■       =  h    >:  >:        hL   N,?n  ,    ,<ik|  |lm> 

ij  it        i  i    i  i      i    km  km         '  ' 

J  k ,  1 ,  m  k    ,  1    ,  m 

r  6i  i  ifin  .5      .  (f,  n  i     i  iVfi  i    ) 

i     kk      II '  mm  v   k  1  m   '      klnr 


(E+e,  -Ei  -e    )  (E+c, -e-,-e    )  (E+e,  ,-e-,  ,  -e 

v       k     1     nr  v       k     1     nr v       k       1 


m" 
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(VivivW(WvW 

k,A,m  T"E+Ek-erem,(E+VEX-S,(E+ek'-ET-em'T 
+   .    .    .  ]   <1V| |jk'>   .  (4.23) 

As  was  done  in  Section  3.3,  the  terms  in  Eq.  (4.23)  can  be  simplified  by 
first  restricting  the  summations  over  all  spin  orbitals  to  summations 

over  occupied  or  unoccupied  spin  orbitals  such  that  the  occupation 

\.      i, 

number  factor  N2,  N,2,,,  ,  is  nonvanishing.  Doing  this  in  the  first  term 
k im  k  I  m  a      ° 

of  Eq.  (4.23)  and  then  summing  over  the  delta  functions  yields 

h     y        <ial|Pq><pq[[Ja>        k       ^  (4  ?/j 

a,p,q      'a  p  q; 


*  * h  try       +  n  <■'■'» 

p,a,b      p  a  b' 


which  are  the  same  second-order  self-energy  diagrams  as  obtained  in  Eqs. 
(3.66)  and  (3.67). 

Restricting  and  spin  orbital  summations  in  the  second  term  of  Eq. 
(4.23)  yields  the  following  expressions 

V      y      y  <ial  IPq>    (f  I  v-F    )  <rsl  lJb>  (a    oa\ 

'a.p.q  b.r.s   ^aV?  (fb,slVfaW>   (B,^)  (4-26) 

+v     y  y  <ial  |pq>       (f       iyf       n       <bc  1  [jr>  ,.   ??, 

\p.,  r.b.c   <^a!Veq»   '^W   ^vW  (4'27) 

+'      r  r  <1"Pl  lab>       (f       I v-f       1       <5rJJj'c> 

2p,a,b  d.q.r   KV?     W    (E+^-e,.)  (4.28) 

+h     T  T        ,    <ipllab>    .     (f        1  v-f        )   ..    <cd|  |  jq>  ,  , 

n  i'   k  „  i  ^    [E+c   -e   -e.  )    ^qcd|VTpabj  TE+e   -c   -eJ  [    '    yj 

p,a,b  q,c,d  p     a     b'       M  p         v       q     c     d' 

Now  substituting   Eq.    (4.18)   for  the  effective   interaction  matrices,  we 
find  that  the  delta   functions   in   Eq.    (4.18)   further  restrict  the  spin 
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orbital  summations  in  such  a  way  that  only  expressions  (4.26)  and  (4. 29' 
are  nonvanishing.  After  some  simplification,  the  nonvanishing  contribu- 
tions are  found  to  be 


h  i  „  a  .  (^l-^fel^)    r\        '4-3° 


a         p,q,r,s  a     p     q' x       a     r     s 


I 


v         <ial  [pqxbgl  [arxprl  |  jb>  ,  . 


a,b         p,q,r     v       a     p     q!K       b     p     r 


E         <1|^<bpllar><rq||Jb>  (432) 

a,b         p,q,r      [t  Ga  ep  eqM        b     r  V 


<ip|lab><balldcxcd||,1p>  r-1T  (        , 

a.b.c.d        P        (E+£p'ea-£b)(E+VVedT 


+1  <ip||ab><qbllpc><ac||jq>  (  } 

a,b,c        p,q      (!Wb}WqVcy 

a,b,c         p,q       U  ep  Ga  VU  eq  Gb  V 

Expressions  (4.30)  and  (4.33)  are  the  only  ones  which  represent  valid 
third-order  diagrams  as  written,  however,  by  interchanging  dummy  indices 
p  +■*■   q  in  Eq.  (4.32)  and  a  ^+  b  in  Eq.  (4.35),  expressions  (4.31)  and 
(4.32)  can  be  combined  to  yield 

'   'I 


I  7, 

a ,b    p,q  ,r 


<ia||pqxbqllarxPrHJb>    L     U 


and  expressions  (4.34)  and  (4.35)  combined  to  yield 

,,■■_„„„■, .     "  "0.. 

a,b,c    p,q 


T     (E^-f-fll^r"1^)      f      A  <4-37> 

,b,c    p,q   (t  Ep  ea  Eb,(E+eq  S  V    _A .17 
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which  now  correspond  to  the  two  third-order  ring  diagrams  as  indicated. 
These  results  verify  that  one  of  our  original  objectives,  which  was  to 
include  all  third-order  ring  and  ladder  diagrams,  has  been  achieved. 

The  diagrammatic  analysis  of  the  third  term  in  Eq.  (4.23)  proceeds 
in  the  same  way  as  that  of  the  first  two  terms,  but  since  the  effective 
interaction  matrix  appears  twice,  it  involves  considerably  more  algebra. 
For  this  reason,  we  simply  display  the  resulting  diagrams  in  Fig.  4. 
It  is  significant  to  realize  that  the  fourth-order  diagrams  in  Fig.  4 
include  not  only  ring  and  ladder  diagrams  but  also  mixed  diagrams  which 
consist  of  both  ring  and  ladder  parts.   In  the  third-order  analysis,  the 
first  term  in  Eq.  (4.18)  was  responsible  for  yielding  the  ladder  diagrams 
while  the  second  and  third  terms  yielded  the  ring  diagrams.   If  we  there- 
fore denote  the  first  term  as  a  ladder  part  and  the  second  and  third 
terms  as  ring  parts,  the  mixed  diagrams  in  fourth  order  dre   found  to 
arise  from  the  product  of  a  ladder  part  and  a  ring  part.   Inducing  the 
results  of  the  fourth-order  analysis  to  higher  orders,  we  conclude  that 
our  second  objective,  which  was  to  sum  all  ring  and  ladder  diagrams  in 
all  orders  of  the  self-energy  expansion,  has  been  exceeded:  not  only 
are  all  ring  and  ladder  diagrams  included  in  the  2p-h  TDA  self-energy, 
but  also  the  mixed  diagrams  which  exhibit  both  ring  and  ladder  parts. 

The  diagonal  2p-h  TDA  self-energy  may  also  be  analyzed  diagrammat- 
ically.  This  analysis  is  even  simpler  than  for  the  full  2p-h  TDA  since 
the  denominator  shifts  are  now  scalars  rather  than  matrices.  A  compar- 
ison of  diagrams  obtained  with  the  denominator  shift  in  Eq.  (4.30) 
versus  that  in  Eq.  (4.22)  will  reveal  the  significance  of  the  factor  of 
'i.  Considering  the  approximation  in  Eqs.  (4.21)  and  (4.22)  first,  we 
obtain  the  following  third-order  expressions 


Figure  4.  Fourth-order  self-energy  diagrams  arising  from 
the  2p-h  TDA. 


100 


'G    TG 
....3 1_.0 

0...  TIL 

..:.:o 

'0..    TO. 


<*-r 


..J> 


101 


!4     i        <ia|  Ipqxpql  |  pqxpq  |  jja> 
4a,P,q  (E+W£q)2 


.j,     v       <ip[  |abxab|  |ab»<ab|  |  jp> 


L       <ia[  lpqxap|  |apxpq|  |ja> 

a,p,q  (E+e   -e   -e    )2 

a      p     q 


,  -^.3q 


(4.38) 


(4.39) 


(4.4o; 


+   E   <ipl  [abxpb[  [pbxabl  |jp> 
a.b.p      (E+VVeb)2 


(4.41) 


The  only  difference  between  these  diagrams  and  the  third-order  diagrams 
of  the  full,  2p-h  TDA  is  that  the  incoming  lines  on  the  middle  interac- 
tion line  have  the  same  labels  as  the  outgoing  lines.  Now  analyzing  the 
approximation  in  Eqs.  (4.19)  and  (4.20),  we  obtain  the  following  third- 
order  expressions 


jj  E   <ia|  |  pqxpq  |  Ipqxpql  [ja> 
2a,P,q      (E+VVV2 


(4.42) 


__,,     y       <ip|  |ab><ab|  |abxab|  | jp> 
2a,b,p  (E^p-^-%)2 


14.43) 


t. 


na  I  I  pax-ap |  |apxpq|  |  ja> 

2 
a,p,q      (E+e  -e  -e  ) 
K  '      v   a  p  q; 


(4.44) 


+   v   <ip|  jabxpbl  [pbxab|  [jp> 
a.b.p      (E+Vca-£b)2 


(4.45; 
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Expressions  (4.44)  and  (4.45)  are  identical  to  (4.40)  and  (4.41)  respec- 
tively; however,  expressions  (4.42)  and  (4.43)  both  differ  by  a  factor 
of  h   from  (4.38)  and  (4.39).  This  discrepancy  is  a  direct  consequence 
of  the  missing  factor  in  the  denominator  shift,  Eq.  (4.20),  and  leads 
to  an  overcounting  of  these  third-order  diagonal  ladder  diagrams  since, 
by  rule  5  in  Table  5,  there  should  be  a  factor  of  H   for  each  pair  of 
equivalent  lines.  Similarly,  it  is  rather  easy  to  show  that  this  approx- 
imation overcounts  all  higher  order  diagrams  containing  this  diagonal 
ladder  part. 

4.4  Computational  Applications  and  Evaluation  of  the  Diagonal  2p-h  TDA 
Self-Energy  ~ "  ' 

The  main  attraction  of  the  diagonal  2p-h  TDA  self-energy  for  the 
calculation  of  ionization  energies  and  electron  affinities  is  its  pseudo 
second-order  structure.  Computational  experience  with  the  diagram  con- 
serving decouplings  has  taught  us  that  the  second-order  self-energy 
approximant  is  both  easily  constructed  and  evaluated.  The  diagonal  2p-h 
TDA  requires  only  the  additional  evaluation  of  Eq.  (4.22)  which  merely 
shifts  the  second-order  self-energy  poles.  As  a  consequence  the  diagonal 
2p-h  TDA  self-energy  mimics  the  exact  self-energy  by  possessing  only 
simple  poles.  Another  consequence  of  the  pseudo  second-order  structure 
is  that  the  energy  dependence  will  have  the  simple  analytic  form  illus- 
trated in  Fig.  2.  This  property,  which  was  absent  in  the  [1,1]  self- 
energy  approximant,  simplifies  the  pole  search  for  the  electron  propaga- 
tor.  In  spite  of  the  simplicity  of  its  pseudo  second-order  structure, 
however,  the  diagonal  2p-h  TDA  self-energy  incorporates  diagonal  ring, 
ladder,  and  mixed  diagrams  in  all  orders  of  the  self-energy  expansion 
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as  was  shown  in  the  previous  section.  This  property  encourages  specula- 
tion that  significantly  more  relaxation  and  correlation  will  be  accounted 
for  in  this  decoupling  than  with  the  second-order  decoupling,  but  the 
accuracy  of  the  corresponding  ionization  energies  or  electron  affinities 
can  only  be  evaluated  via  actual  computational  applications.  As  was 
done  with  the  diagram  conserving  decouplings,  ionization  energies  for 
the  water  and  formaldehyde  molecules  were  computed  using  the  diagonal 
2p-h  TDA  self-energy. 

For  the  water  molecule,  calculations  were  performed  with  both  the 
14  and  26  contracted  Gaussian  orbital  basis  sets  described  in  Chapters 
2  and  3.  The  principal  ionization  energies  computed  with  the  diagonal 

2p-h  TDA  and  the  diagonal  2p-h  TDA  plus  third-order  constant  energy 

(2)         (2)        (3) 
diagrams  (denoted  e(E)shjft  and  Z(E)SHIFT  +  Z [<*>))   are  presented  in  Table 

11.  Comparing  these  results  with  those  in  Table  9  for  the  second-order 
and  [1,1]  self-energy  approximants  reveals  that  each  ionization  has  been 
shifted  to  higher  energy.  This  shift  has  led  to  a  significant  improve- 
ment in  the  valence  ionization  energies  (3a,,  lb.,  and  lb?)  which  are 
now  within  approximately  0.5  eV  of  the  experimental  results.  For  the 
inner  valence  (2a,)  and  core  (la.)  ionizations,  however,  this  energy 
shift  leads  to  worse  agreement.  In  addition,  the  diagonal  2p-h  TDA 
results  for  the  la,  ionization  exhibit  an  enormous  basis  dependence. 
The  addition  of  polarization  functions  in  the  26  orbital  basis  has 
yielded  nearly  13  eV  in  additional  relaxation.  The  most  probable  explana- 
tion for  this  basis  dependence  is  that  the  2p-h  TDA  self-energy  poles 
are   determined  by  Hartree-Fock  orbital  energies  and  two-electron  inte- 
grals rather  than  by  orbital  energies  alone  as  with  the  second-order 
self-energy.  The  orbital  energies  are  rather  insentitive  to  basis 
changes,  whereas  the  two-electron  integrals  are  not. 
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Calculations  for  the  formaldehyde  molecule  were  performed  with  the 
24  and  42  orbital  basis  sets  described  in  Chapter  3.  Ionization  energies 
were  computed  using  the  diagonal  2p-h  TDA  and  are  presented  in  Table  12. 
The  third-order  constant  energy  diagrams  were  not  evaluated  for  this 
molecule.  Similar  to  the  water  results,  the  diagonal  2p-h  TDA  results 
for  formaldehyde  are  also  consistently  higher  in  energy  than  the  second- 
order  results  (Table  9).  This  shift  considerably  improves  the  valence 
ionization  energies;  however,  the  average  deviation  from  the  experimental 
results  remains  approximately  0.8  eV.  The  core  ionizations  within  the 
diagonal  2p-h  TDA  suffer  a  small  deterioration  in  accuracy  but  do  not 
exhibit  the  extreme  basis  dependence  which  was  observed  in  the  water 
calculations.  Part  of  the  discrepancies  between  the  diagonal  2p-h  TDA 
and  the  experimental  results  can  certainly  be  eliminated  by  further  basis 
saturation;  however,  in  the  next  chapter,  we  will  propose  that  even  ioni- 
zation energies  of  1.0  eV  accuracy  are  usually  sufficient  to  unambiguously 
interpret  photoelectron  spectra  if  combined  with  a  calculation  of  rela- 
tive photoionization  intensities. 

Cederbaum  and  co-workers  have  recently  developed  computer  programs 
which  implement  the  full,  nondiagonal  2p-h  TDA  to  the  self-energy  and 
have  reported  several  molecular  applications  (Cederbaum  et  al_. ,  1977, 
Schirmer  et  aj_. ,  1977,  Cederbaum  ejt  aj_. ,  1978,  and  Schirmer  et  al . , 
1978).   In  these  calculations,  they  claim  only  a  1.0  eV  accuracy  and  rely 
heavily  on  vibrational  analyses  to  assist  with  the  interpretation  of 
photoelectron  spectra.  The  off-diagonal  matrix  elements  seem  to  have 
little  importance  in  the  valence  region  because  the  propagator  poles  are 
relatively  well-separated.  In  the  inner  valence  and  core  regions,  how- 
ever, where  principal  ionization  poles  and  shake-up  poles  overlap  and 
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Table  12.  Formaldehyde  Results  Obtained  with  the  Diagonal  2p-h 
TDA  Self-Energy. 


24CGT0's  42  CGTO's 


(2)  (2) 

Koopmans   i(EL,,T[--r   Koopmans   E(E),-utct     Exp. 


lal 

560.12 

542.86 

559.81 

542.09 

539. 43a 

2a, 

309.09 

299.82 

308.87 

299.31 

294. 21a 

3a, 

38.94 

35.60 

38.18 

34.31 

34. 2b 

4a  j 

23.39 

21.75 

23.30 

21.79 

21.15b 

5aj 

17.29 

14.99 

17.38 

15.33 

16. 2C 

lb. 

14.56 

14.11 

14.45 

14.03 

14.5C'd 

lb2 

19.47 

17.91 

19.08 

17.77 

17. 0C 

2b0 

12.06 

9.92 

11.93 

10.11 

10. 9C 

aJolly  and  Schaaf  (1976). 

bHood  et  aj_.  (1976). 

Estimated  center  of  gravity  (Cederbaum  and  Domcke,  1977)  from 
spectrum  of  Turner  et  a]_.  (1970). 

d14.38(8)  VIP  (Turner  et  al.,  1970). 
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interact,  level  shifts  and  intensity  changes  are  observed.   In  these 
regions,  even  the  nondiagonal  2p-h  TDA  is  not  fully  satisfactory  since 
ion-state  relaxation  and  hole-hole  interactions,  neither  of  which  are 
described  by  this  self-energy  approximation,  may  also  be  important 
(Wendin,  1979). 


CHAPTER  5 
PHOTOIONIZATION  INTENSITIES 


5.1  Introduction 


The  evaluation  of  each  decoupling  approximation  in  the  preceding 
chapters  was  based  on  the  comparison  of  propagator  poles  to  experimental 
ionization  energies.  This  criterion  represents  a  particular  bias  since 
it  does  not  reflect  the  quality  of  the  Feynman-Dyson  amplitudes  (defined 
in  Section  1.1).  The  Feynman-Dyson  amplitudes  determine  the  spectral 
density  function  (Linderberg  and  Ohrn,  1973) 

£  fk(x)f*(x')S(E-Ek)        E  >  u 
k 
A(x,x ' ;E)  = 

'I   9k(x)g"(x')6(E-Ek)        E  <  u         (5.1) 

which  contains  a  plethora  of  useful  information.  This  is  evidenced  by 
the  relation  of  the  spectral  density  to  the  single-particle,  reduced 
density  matrix  (Linderberg  and  Ohrn,  1973) 


Y(x.x')  =    A(x,x';E)dE  .  (5.2) 

It  is  important,  therefore,  to  choose  a  decoupling  approximation  which 
not  only  yields  accurate  ionization  energies  but  also  an  accurate  spectral 
density.  The  quality  of  the  spectral  density  is  somewhat  more  difficult 
to  evaluate  since  there  are  no  experimental  data  with  which  it  can  be 
directly  compared.  One  evaluation  procedure,  however,  might  involve  the 
calculation  of  averages  of  specific  one-electron  operators  from  the 
reduced  density  matrix.  The  averages  could  then  be  compared  with 
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experimental  results.  Another  procedure  which  is  more  closely  related 
to  our  goal  of  interpreting  photoelectron  spectra  is  the  calculation  of 
photoionization  intensities  or  cross-sections.  A  theoretical  prediction 
of  relative  photoionization  intensities  can  simplify  orbital  assignments 
when  the  ionization  energies  are  not  accurately  given.  A  theoretically 
predicted  variation  in  relative  intensities  with  photon  energy  is  par- 
ticularly useful  if  photoelectron  spectra  are  available  with  different 
ionization  sources  (Katrib  et  aj_.  ,  1973).   In  the  following  sections, 
we  derive  computational  expressions  which  relate  the  Feynman-Dyson  ampli- 
tudes to  the  total  photoionization  cross-section,  discuss  the  major 
approximations  assumed,  and  then  present  several  applications. 

5.2  Derivation  of  Computational  Formulae  for  the  Total  Photoionization 
Cross-Section 

The  differential  cross-section  for  photoionization  derived  from 

first-order,  time-dependent  perturbation  theory  using  a  semi-classical 

model  for  the  interaction  of  radiation  and  matter  is  (Bethe  and  Salpeter, 

1957,  Kaplan  and  Markin,  1968,  Smith,  1971) 


da         4tt    Ik  J 

s_  _  '    f ' 

f       c|Ap.|   co 


<N,s|Z  t.   •  $.     |N> 
k     K       K 


(5.3) 


In  this  equation,  kf  is  the  momentum  of  the  ejected  photoelectron  and 

A  =  I   A,  is  the  vector  potential.  For  a  closed-shell  system,  the  initial 

k  K 
state  |N>  can  be  represented  by  an  anti symmetrized  N-electron  wavefunction 

|N>  =  $0(x1,x2,  .  .  .  xN)  (5.4) 

and  the  final  state  |N,s>  is  represented  by 
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|N,s>  =  (N/2)'5  0AS  [v(tf,r)a(c)*se(x1,x2,  .  .  .  x^j) 

-  v(tf)f)3(?)$sa(xrx2,  .  .  .  x^)].     (5.5) 

Here  $   and  $   are  the  two,  doublet  spin  components  of  the  (N-l)- 
electron  ion,  v(lL,r)  denotes  the  wavefunction  of  the  photoelectron,  and 
Or,,-   is  an  antisymmetrizer 

°AS  =  N~1^   \lt   PkN]  •  (5-6) 

The  form  of  Eq.  (5.5)  guarantees  that  the  singlet  spin  symmetry  of  the 
system  is  preserved. 

Evaluating  the  matrix  element  in  Eq.  (5.3)  yields  (Purvis  and  Ohrn, 
1975a) 


<N. 

k 


1-sfE  Ak-  Vk  |N>  =  /I  |  v*(tffr)^s-  ^  gs(r)d^ 

+  /?  J  V*(tf,f)ps(f)dr  (5.7) 

where 

N'%gs(r)x(c)  =  |  0*  (xr  .  .  .  xN_1)$0(xr  .  .  .  xN_i;r,x(c)) 

x  dXj  .  .  .  dxN_}  (5.8) 

and 

N_isPs(r)x(c)  =  (N-l)  |  $*  (xp  .  .  .  Y1)^-V1$0(x1,  .  .  .  xN_r- 

r,x(c))dXj  .  .  .  dxN_1  .  (5.9) 

The  first  term  in  Eq.  (5.7)  relates  the  Feynman-Dyson  amplitude  g  (r)  to 
the  photoionization  cross-section  and  the  second  term  arises  from  the 
nonorthogonal  ity  of  v(kf.,r)  and  $«.  When  v(£f,~r)  is  strongly  orthogonal 
to  $0,  this  term  vanishes.  Even  when  v(k^,f)  is  not  strongly  orthogonal 
to  $0,  the  first  term  in  Eq.  (5.7)  will  dominate  if  the  kinetic  energy 
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of  the  photoelectron  is  much  greater  than  its  binding  energy  (Rabalais 
et  aj_. ,  1974),  therefore  we  shall  neglect  this  term: 


da 
s 


3ir2|t 


c|Ar 


v*(tf,r)Xs-^sgs(r)df 


5.10! 


Further  simplification  can  be  obtained  by  neglecting  retardation  of 
the  photoelectron  momentum  by  the  photon  momentum  (also  called  the  dipole 
approximation  in  analogy  with  photon  induced  transitions  in  the  discrete 
spectrum,  Bethe  and  Salpeter,  1957,  Steinfeld,  1974).  The  vector  poten- 
tial As  has  the  following  plane  wave  decomposition  in  terms  of  the  inci- 
dent photon  momentum  k  and  polarization  n 

/ts  =  \tQ\   exp  (il<u  •  rs)n  .  (5. 11) 

If  the  wavelength  of  the  incident  photon  is  large  compared  with  molecular 
dimensions,  the  exponential  may  be  approximated  by  the  first  term  in  its 
Taylor  series  expansion  (unity) 

A 


|A0|n 


5.12 


With  this  approximation,  the  differential  photoionization  cross-section 
becomes 


das 

d£L 


where 


?  = 


Btt2  I  ]<, 


n 


v   (kf,r)Vgs(r)dr   . 


(5.13) 


(5.14) 
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Owing  to  the  random  orientation  of  molecules  in  a  gaseous  sample, 
the  experimentally  observed  photoionization  intensities  in  the  molecular 
reference  frame  represent  an  average  over  all  incident  photon  directions. 
Furthermore,  if  the  incident  photon  beam  is  unpolarized,  we  must  also 
average  over  photon  polarizations.  Making  the  appropriate  averages  in 
Eq.  (5.13)  (Smith,  1971),  we  obtain 


do 


[  8tt2||<  |^      , 


Since  the  two  polarization  directions,  n,  and  rL,  and  the  incident  photon 
direction,  ky[kj,  form  a  right-handed  system  of  axes,  we  can  write 

i^i2  =  iv?i2  +  iv^i2  +  iV^i2  /  ^f  ■  <5-16) 

and  Eq.  (5.15)  becomes 


[-nr-J  mj^2-\K-f\2ntjh^        (5.i7) 


da 
s_ 


(1  -  cos2  9  }dfi  (5.18) 

WO)  v     ' 

da      8tt2|  ft  f 

d^  =   -EST1"  HT  •  (5.19) 

In  order  to  evaluate  |f*|  ,  some  form  for  the  photoelectron  wave- 
function  v("kf,r)  must  now  be  chosen.   In  principle,  the  photoelectron 
wavefunction  could  be  obtained  by  the  solution  of  the  Bethe-Salpeter 
equation  for  the  polarization  propagator  where  the  superoperator  resol- 
vent has  been  modified  to  include  the  time-dependent  interaction  of  the 
radiation  and  matter  fields  (see  e.g.  Csanak  et  aJL  ,  1971).  A  solution 
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of  this  type  would  require  the  use  of  continuum  functions  as  well  as 
discrete  functions  in  the  molecular  basis  and  is  not  yet  feasible  owing 
to  several  formal  and  practical  difficulties.  Alternatively,  we  seek  a 
simple  but  accurate,  analytic  representation  of  the  photoelectron  wave- 
function.  For  photoionization  of  atoms  or  molecules  with  high  (tetra- 
hedral  or  octahedral)  symmetry,  the  electronic  potential  of  the  ion  is 
nearly  spherically  symmetric,  and  v(tf,f)  may  be  asymptotically  repre- 
sented by  a  plane  wave  plus  incoming  Coulomb  waves  (see  e.g.  Smith,  1971). 
For  molecules  with  lower  symmetry,  distortions  of  the  electronic  poten- 
tial enormously  complicate  the  nature  of  the  incoming  waves.   In  this 
derivation,  the  incoming  waves  are  neglected,  and  the  photoelectron  is 
simply  represented  by  the  plane  wave  part 

v(l<f,r)  =  (2tt)"3/2  exp{itf  ■   r)  .  (5.20) 

The  applicability  and  implications  of  this  approximation  will  be  dis- 
cussed in  the  following  section.  With  this  choice,  Eq.  (5.14)  can  be 
integrated  by  parts  and  yields 

?  =  (2^)"3/2  |  gs(r)v^exp(-i£f  •  r)dr  (5.21) 

=  -i£f(2^)"3/2  I  exp(-Htf  •  r)gs(r)dr  .  (5.22) 

In  our  computational  scheme,  the  Feynman-Dyson  amplitudes  in  Eq. 
(5.22)  are  represented  by  a  linear  combination  of  Hartree-Fock  orbitals. 
The  Hartree-Fock  orbitals  can  be  decomposed  into  contracted  Cartesian 
Gaussian  functions  on  each  atomic  center,  and  each  contracted  Gaussian 
function  can  be  further  decomposed  into  a  sum  of  primitive  Gaussian 
functions.  Ultimately  therefore,  the  Feynman-Dyson  amplitudes  can  be 
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represented  as  some  linear  combination  of  primitive  Gaussian  functions 
on  each  atomic  center 


gs(r)  =  I     c^  ^^  (r  -  ft  J  .  (5.23! 

a,k 


'ak  ak  •     ,</ 


Here  a  represents  the  sum  over  atomic  centers  and  k  the  sum  over  primitive 
Gaussian  functions  on  each  center.  Transforming  the  photoelectron  posi- 
tion vector  to  the  coordinate  frames  of  each  atomic  center  by  the  sub- 
stitution r'  =  r-R  ,  Eq.  (5.22)  can  be  re-expressed  in  terms  of  the  prim- 
itive computational  basis 


ft  =  -itf  (2tt)"3/2  I  exp(-itf  •  ft  )c^k  j  exp(-itf  •  r1  J^fr'  )dr'  . 

(5.24) 
The  integral  over  r'  in  Eq.  (5.24)  represents  a  Fourier  transform  of 
each  primitive  Gaussian  function,  and  formulae  for  its  evaluation  are 
derived  by  Kaijser  and  Smith  (1977). 

The  final  step  in  this  derivation  is  to  integrate  the  differential 
photoionization  cross-section  (Eq.  (5.19))  over  all  photoelectron  direc- 
tions in  the  solid  angle  dfif 

8tt  kf  r         ? 

Since  P  involves  a  sum  over  atomic  centers,  Eq.  (5.25)  will  contain  both 
one-  and  two-center  contributions  (Schweig  and  Thiel,  1974) 

as  =  4c7T  \       *     Cak  C61  <&,B1  (5-25^ 

a,k  6,1 

The  one-center  terms  (a=B)  have  the  form 
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>s 


Q \     t   =       $  u(M$  i(k*)dn#  (5.27) 

ak,al        J      ak     f     al      f       f  w,t'' 

where  $  .  (k^)   represents  the  Fourier  transform 

*ak(tf)   »   (2tt)_3/2  J  exp(-itf  •  r'  )*ak(r'  )dr '  (5.28) 

and  the  two-center  terms   (a^3)  are  given  by 

Qak,61    =  j<k<V  VVe*P^V  U^f  (5-29^ 

where  R   =  R  -  R    Owing  to  the  orthogonality  of  the  spherical  har- 
monics which  describe  the  angular  dependence  of  the  Fourier  transforms, 
the  one-center  terms  are  easily  evaluated.  The  two-center  terms  are 
complicated  slightly  by  the  exponential  factor,  but  these  too  can  be 
evaluated  analytically.  For  this  purpose,  it  is  convenient  to  define 
the  integral 


I(1iVW  =  j  ^VM%m2(V¥exP^V  VdV   (5-30) 

where  /-j  m  and  V-\?m?   are  real  spherical  harmonics  (Harris,  1973)  repre- 
senting  the  angular  dependence  of  the  transforms  $  ,  and  $_,  respectively. 
Using  the  expansion 

exp(±ltf  •  ftaB)  -  4^Q  ^^^VMaB^m^VV9^'  ^31> 

Eq.  (5.30)  becomes 

m        1     -,         m..m„-m 
1(1  m  l?m  )  =  4tt  l        l    (-i )  'j  (kfR  JC,  V  ,  ylm(e,rf>).     (5.32) 
1         1=0  m=-l  ''   '12 

In  these  equations,  {$*,§„)   represents  the  photoelectron  direction  in 

the  atomic  reference  frame,  ( B , cb )  represents  the  direction  of  R  „in  the 

r  aB 

molecular  frame,  j -.  ( k^R  )  is  a  spherical  Bessel  function  of  1-th  order, 

m]m2-m 
and  C]  i  i  is  a  Clebsh-Gordan  coefficient  defined  by  (Harris,  1973) 
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m,m„m 


\izi  "  J  \m^r^yi^r^yiJ^f)^f   ■  (5.33) 

When  the  distance  vector  ft   coincides  with  the  molecular  z-axis,  the 
expansion  in  Eq.  (5.31)  simplifies  to 

exp(±ltf  -fta0)  =  (4-rr)-  F,    (±1 )]  (21+1)^  (y^g)^^  ,*f )  (5.34) 


1 

and  Eq.  (5.30)  becomes 


1 1  -,  +  l^j 

i,       1   ^     i     i  m1m90 

1(1lml',2",2»-<*'>\m2  ,.„*_ ,  i'""  (21+1)SCl}l2l  VkfV 

(5.35) 
For  arbitrary  directions  of  ft   which  do  not  coincide  with  the  molecular 
z-axis,  a  new  coordinate  frame  can  be  defined  in  which  ft  „  does  coincide 
with  this  axis,  and  the  spherical  harmonics  can  be  transformed  to  this 
coordinate  frame  using  the  familiar  rotation  matrices  (Schweig  and  Thiel , 
1974).   In  this  way,  Eq.  (5.35)  can  be  used  to  evaluate  all  two-center 
integrals  arising  from  any  molecular  geometry. 

5.3  Discussion  of  Approximations 

Relatively  little  work  has  been  devoted  to  the  theoretical  calcula- 
tion of  photoionization  cross-sections  for  molecules  compared  with  that 
for  atoms  (see  e.g.  Marr,  1967,  Steward,  1967,  Kelly,  1976).  The  major 
impediments  until  recently  have  been  a  lack  of  sufficiently  accurate 
molecular  wavefunctions  and  an  absence  of  accurate  analytic  representa- 
tions for  the  photoelectron  wavefunction  (Kaplan  and  Markin,  1967, 
Schweig  and  Thiel,  1974).  With  the  development  of  efficient,  molecular 
integral  and  Hartree-Fock  programs,  Hartree-Fock  wavefunctions  are  now 
readily  available  for  a  large  number  of  molecules.  This  availability 
in  turn  has  stimulated  several  theoretical  calculations  of  molecular 
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photoionization  cross-sections  using  the  frozen  orbital  approximation 
(Rabalais  et  aj_. ,  1974,  Dewar  et  al_. ,  1975,  Allison  and  Cavell,  1978, 
and  Cavell  and  Allison,  1978). 

In  the  frozen  orbital  approximation,  the  N-electron  reference  state 
is  assumed  to  be  the  Hartree-Fock  qround  state,  and  the  ion  states  are 
constructed  by  removing  the  orbital  corresponding  to  the  ionized  elec- 
tron. This  approximation  neglects  both  correlation  corrections  in  the 
ground  and  ion  states  and  ion  state  relaxation  (see  Introduction). 
Following  Purvis  and  Ohrn  (1975a),  we  have  replaced  the  frozen  orbital 
approximation  by  a  many-electron  treatment  which  incorporates  both  relax- 
ation and  correlation  corrections.  This  treatment  derives  from  the  use 
of  the  Feynman-Dyson  amplitudes  obtained  from  the  electron  propagator  to 
compute  the  photoionization  cross-section. 

The  form  of  the  photoelectron  wavefunction  is  still  a  major  problem 
in  the  calculation  of  molecular  photoionization  cross-sections  and  repre- 
sents the  most  critical  step  in  our  derivation.  The  plane  wave  approxi- 
mation was  chosen  for  its  simplicity,  however,  it  has  several  serious 
limitations.  The  most  serious  limitation  is  its  failure  to  correctly 
predict  experimentally  observed  angular  distributions  (Bethe  and  Salpeter, 
1957,  Schweig  and  Thiel,  1974).  This  deficiency  is  not  readily  apparent 
when  the  differential  cross-section  is  averaged  over  all  photoelectron 
directions,  and  in  our  computational  applications,  we  compute  only  spher- 
ically averaged,  total  cross-sections.  Lohr  (1972)  has  shown  that  by 
retaining  the  second  term  in  Eq.  (5.7),  qualitatively  correct  angular 
distributions  may  be  obtained.  The  retention  of  this  term  is  equivalent 
to  orthogonalizing  the  plane  wave  to  every  one-electron  function  in  the 
N-electron  ground  state  and  is  known  as  the  orthogonal ized  plane  wave 
(OPW)  approximation. 
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A  second  limitation  of  the  plane  wave  approximation  is  the  implicit 
neglect  of  electrostatic  interactions  between  the  photoelectron  and  the 
molecular  ion.  For  ionization  processes  near  threshold  where  the  photo- 
electron  leaves  with  low  kinetic  energy,  these  interactions  are  espe- 
cially important,  and  the  wavefunction  of  the  photoelectron  exhibits  a 
decrease  in  wavelength  as  r  ->  fj.  The  OPW  approximation  again  partially 
corrects  this  deficiency  (Lohr,  1972)  but  exhibits  a  rather  abrupt  change 
in  wavelength  which  is  more  characteristic  of  a  short-range  potential 
rather  than  the  long-range  Coulomb  potential.  Other  representations  of 
low  energy  photoelectrons,  e.g.  multicentered  Coulomb  wave  expansions, 
have  also  been  proposed  (Tuckwell,  1970).  As  the  kinetic  energy  of  the 
photoelectron  tends  to  higher  energies,  the  effect  of  electrostatic 
interactions  on  the  total  cross-section  becomes  negligible,  and  the  plane 
wave  approximation  becomes  sufficiently  accurate  (Bethe  and  Salpeter, 
1957). 

Unfortunately,  as  the  plane  wave  approximation  becomes  more  suitable 
at  high  photoelectron  energies,  the  dipole  approximation  (Eq.  (5.12)) 
deteriorates  and  must  be  re-examined.  Bethe  and  Salpeter  (1957)  have 
shown  that  when  retardation  effects  are  included  in  the  calculation  of 
differential  photoionization  cross-sections,  the  lowest  order  correction 
is  proportional  to  (|v|/c),  where  v  is  the  photoelectron  velocity.  This 
correction,  however,  only  effects  the  angular  distribution  and  vanishes 
when  the  differential  cross-section  is  averaged  over  all  incident  photon 
directions.  Higher-order  corrections  for  retardation  are  proportional 
to  (|v|/c)  as  are  the  relativistic  corrections,  but  these  are  usually 
negligible  even  when  the  photon  wavelength  is  comparable  to  the  molecular 
dimensions. 
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5 . 4  Computational  Applications 

In  this  section,  the  relative  photoionization  intensities  for  the 
water  and  acetylene  molecules  are  computed  using  the  Feynman-Dyson  ampli- 
tudes obtained  in  the  second-order,  diagram  conserving  and  the  diagonal 
2p-h  TDA,  renormalized  decoupling  approximations.  Comparative  calcula- 
tions corresponding  to  a  Mg  K  photon  source  (tioj  =  1253.6  eV),  for  which 
the  plane  wave  approximation  is  expected  to  be  good,  and  a  He  (II)  photon 
source  (nw  =  40.81  eV),  for  which  the  plane  wave  approximation  is  expected 
to  be  less  accurate,  are   presented.  These  results  are  further  compared 
with  intensities  obtained  using  the  frozen  orbital  approximation  in  order 
to  assess  the  magnitude  of  many-electron  correlation  and  relaxation  cor- 
rections. For  acetylene,  the  dependence  of  intensity  on  photon  energy 
for  the  valence  orbitals  is  plotted,  and  orbital  and  density  difference 
plots  are  presented  and  discussed. 

The  relative  photoionization  intensities  for  water  corresponding  to 
a  Mg  «a  photon  source  and  a  He  (II)  photon  source  are  presented  in  Tables 
13  and  14,  respectively.  These  results  were  obtained  using  the  26  con- 
tracted Gaussian  orbital  basis  at  the  experimentally  determined  equilib- 
rium geometry  as  described  in  detail  in  Section  3.5  and  are  scaled  rela- 
tive to  the  3aj  intensity.  As  expected,  the  relative  intensities  com- 
puted for  the  Mg  K  source  in  Table  13  compare  reasonably  well  with  those 
obtained  experimentally  (Rabalais  et  al_. ,  1974).  The  orthogonal ized 
plane  wave  results  of  Rabalais  et  al_.  (1975)  are  presented  for  compari- 
son and  also  correctly  order  the  relative  intensities.   It  is  interesting 
to  note  the  larger  discrepency  between  the  lb,  and  lb?  intensities  in 
our  calculations  compared  to  Rabalais  et  a]_.  than  in  the  2a,  and  3a-, 
intensities.  Since  for  X-ray  photons  the  OPW  corrections  are  not  expected 
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Table  13.  Relative  Photoionization  Intensities  for  Water 
Excited  by  Mg  K  (1253.6  eV). 


Orbital 

F0a 
4.69 

(2) 
2(E) 

1.32 
2.57 

(2) 

£(E); 

SHIFT  + 

(3) 

s(-) 

0PWb 

4.50 

Exp. 

2a , 

4.15 

3.84 

3a-, 

1.00 

1.00 

1.00 

1.00 

1.00 

lbl 

0.62 

0.52 

0.55 

0.18 

0.36 

lb2 

0.40 

0.35 

0.37 

0.09 

0.31 

Frozen  orbital  approximation, 
bRabalais  et  al.  (1974). 
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Table  14.  Relative  Photoionization  Intensities  for  Water 
Excited  by  He(II)  (40.81  eV). 


a    (2)    (2)        (3)     h 

Orbital    F0d   l(E)   z(E)SHIFT  +  E(»)   0PWD      Exp. 


2al 

0.55 

0.28 
0.66 

3a, 

1.00 

1.00 

lbl 

1.16 

1.15 

lb2 

0.83 

0.84 

0.84  0.26  not  observed 

1.00  1.00      1.00 

1.15  0.86      0.96 

0.83  1.10      0.80 


Frozen  orbital  approximation. 
bRabalais  et  al_.  (1975). 
cRabalais  et  al.  (1974). 
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to  be  large,  this  discrepency  is  most  likely  a  basis  set  effect.  Our 
basis  contained  d-type  polarization  functions  on  oxygen  as  well  as  p-type 
polarization  functions  on  the  hydrogen  atoms,  whereas  Rabalais  et  al_. 
have  used  only  a  minimal  basis. 

The  relative  intensities  presented  in  Table  14  do  not  correctly 
match  the  experimentally  observed  intensities.  This  is  not  particularly 
surprising  owing  to  the  inadequacy  of  the  plane  wave  approximation  at 
low  photon  energies.  It  is  more  surprising  that  the  OPW  results  also 
yield  an  incorrect  ordering  of  the  intensities.  The  failure  of  the  OPW 
approximation  indicates  the  necessity  for  a  more  accurate  photoelectron 
wavef unction. 

Many-electron  correlation  and  relaxation  corrections  appear  negli- 
gible in  the  intensities  of  both  Tables  13  and  14.  The  deviations  be- 
tween the  frozen  orbital  approximation  (FO),  the  second-order  self-energy 

(2) 
approximation  (z(E)),  and  the  diagonal  2p-h  TDA  with  third-order  constant 

(2)        (3) 
energy  diagrams  (x(E)SHIFT  +  E(°°)),  are  small  and  may  be  attributed  to 

differences  in  the  photoelectron  momenta.   Different  ionization  energies 
obtained  with  different  decoupling  approximations  were  used  to  compute 
the  photoelectron  momenta  from  the  conservation  of  total  energy 

|kf|  =  [2(w  -  I.P.)]?S  (in  atomic  units)  ,  (5.36) 

consequently,  deviations  in  the  I.P.'s  result  in  deviations  in  |k\-|. 

The  relative  intensities  computed  with  the  diagonal  2p-h  TDA  plus 
third-order  constant  energy  diagrams  for  the  Mg  K  photon  source  are 
plotted  in  Fig.  5.   In  order  to  compare  the  theoretical  spectrum  to  the 
experimental  ESCA  (Electron  Spectroscopy  for  Chemical  Analysis)  spectrum 
(Siegbahn  e_t  al_. ,  1969)  which  is  sketched  as  an  insert,  the  ionization 


Figure  5.  A  plot  of  the  theoretical  ESCA  spectrum  for  the  valence 
ionizations  of  water.  The  experimental  spectrum  of 
Siegbahn  et  a]_.  (1969)  is  sketched  in  the  insert. 


Kii.uft.Lii'i^ii!)  mm  m>mm  \m  mm 


124 


JO  Ti  20  \i 

■.;--h  :-'!'")i.  %> "^   -  f'rcioi         0.'/H!;/.'h         0LUQ. 


125 

lines  were  given  a  Lorentzian  width  estimated  from  the  experimental 
spectrum.  The  major  features  of  the  experimental  spectrum  are  well- 
reproduced  with  the  exception  of  the  weak  peak  at  approximately  23  eV. 
This  peak  has  been  attributed  to  the  ionization  of  2a,  electrons  by  the 
Mg  Kqj    satellite  in  the  photon  source  (Siegbahn  et  al_. ,  1969)  and 
therefore  should  not  appear  in  the  theoretical  spectrum  since  it  was 
calculated  for  a  purely  monochromatic  source. 

A  recent  experimental  study  of  relative  photoionization  intensities 
of  acetylene  with  various  photon  sources  (Cavell  and  Allison,  1978) 
motivated  our  examination  of  this  molecule.  As  preparatory  steps  to  the 
calculation  of  relative  intensities,  Hartree-Fock  and  electron  propagator 
calculations  were  performed  at  the  experimental  equilibrium  geometry, 
R(C-C)  =  2.279  a.u.  and  R(C-H)  =  2.005  a.u.  (Buenker  et  al_. ,  1967),  with 
two  different  basis  sets.  The  first  basis  consisted  of  Huzinaga's  9s, 5p 
primitive  basis  for  carbon  and  4s  basis  (unsealed)  for  hydrogen  (Huzinaga, 
1965)  contracted  with  Dunning's  coefficients  (Dunning,  1970)  to  4s, 2p 
on  carbon  and  2s  on  hydrogen  (see  Table  8).  The  complete  molecular  basis 
consisted  of  24  contracted  Gaussian  orbitals  and  yielded  a  Hartree-Fock 
total  energy  of  E(HF)  =  -76.7948  H.  The  second  basis  augmented  the  first 
by  the  addition  of  a  set  of  d  functions  on  each  carbon  atom  (ex  ,  =  0.60) 
and  a  set  of  p  functions  on  each  hydrogen  atom  (a  =  0.75).  The  addition 
of  these   diffuse,  polarization  functions  brought  the  size  of  the  basis 
to  42  contracted  Gaussian  orbitals  and  yielded  a  Hartree-Fock  total 
energy  of  E(HF)  =  -76.8267  H. 

Valence  ionization  energies  were  computed  with  the  second-order 
self-energy  approximation  and  the  diagonal  2p-h  TDA  for  each  basis.  The 
results  for  the  24  orbital  basis  are   presented  in  Table  15  and  the 
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Table  15.  Valence  Ionization  Energies  for  Acetylene 
(24  CGTO's). 


Orbital 

Koopmans 

(2) 

(2) 
E^SHIFT 

Ced.a 
28.9 

Exp. 

2ag 

- 

38.29 

27.6 

Shake-up 

2o, 

28.42 

25.47 

26.48 

23.9 

23.5 

\ 

18.55 

16.39 

17.08 

16.4 

16.8 

hu 

11.32 

11.20 

11.29 

10.8 

11.4 

2\ 

20.81 

18.21 

19.12 

18.0 

18.7 

aCederbaum  et  al_.  (1978). 
bCavell  and  Allison  (1978) 
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results  for  the  42  orbital  basis  appear  in  Table  16.  In  each  table  the 
full,  nondiagonal  2p-h  TDA  results  of  Cederbaum  et  aj_.  (1978)  and  the 
experimental  results  of  Cavell  and  Allison  (1978)  are  included  for  com- 
parison. In  contrast  to  previous  calculations  reported  here,  the  ioni- 
zation energies  obtained  for  acetylene  are  larger  (with  the  exception  of 
the  1itu  ionization)  than  the  experimental  values.  A  breakdown  in  the 
quasi-particle  picture  for  the  2c  ionization  which  is  evidenced  by  the 
shake-up  pole  may  account  for  the  particularly  poor  results  for  these 
ionizations.  The  off-diagonal  2p-h  TDA  contributions  considerably  im- 
prove these  two  ionization  energies  (Cederbaum  et  aj_. ,  1978),  however, 
the  shake-up  energy  still  disagrees  with  the  experimental  value  by  more 
than  an  electron  volt. 

Relative  photoionization  intensities  for  acetylene  are  presented  in 
Table  17  for  a  Mg  K  photon  source  and  in  Table  18  for  a  He  (II)  photon 

source.   In  both  tables  the  intensities  were  scaled  relative  to  the  Itt 

u 

intensity.  As  for  water,  the  intensities  corresponding  to  the  Mg  K 

source  are  in  reasonable  agreement  with  experimental  intensities  and  show 

only  a  slight  diminution  when  correlation  and  relaxation  effects  are 

included.  The  OPW  result  of  Cavell  and  Allison  (1978)  for  the  2a 

9 

ionization  in  Table  17  seems  unexplicably  high,  and  the  3a  ionization 
is  not  observed  experimentally  (Cavell  and  Allison,  1978). 

The  relative  intensities  for  acetylene  computed  with  the  He  (II) 
source  exhibit  better  agreement  with  the  experimental  results  than  did 
the  water  results.  Here,  only  the  relative  intensities  of  the  two 
weakest  ionizations,  the  2a  and  2a  ,  were  reversed.  The  OPW  results 
predict  the  correct  ordering  but  attribute  a  much  weaker  intensity  to 
the  2a  ionization  than  is  observed  experimentally. 
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Table  16.  Valence  Ionization  Energies  for  Acetylene 
(42  CGTO's). 


Orbital 

Koopmans 

(2) 

(2) 

IE) 
K    ;SHIFT 

34.45 

Ced.a 
28.9 

Exp.b 

2a 
q 

38.20 

27.6 

Shake-up 

\ 

28.07 

24.56 

25.73 

23.9 

23.5 

3a9 

18.46 

16.56 

17.24 

16.4 

16.8 

K 

11.16 

11.13 

11.24 

10.8 

11.4 

2°„ 

20.90 

no  results 

19.47 

18.0 

18.7 

aCederbaum  et  al_.  (1978). 
Cavell  and  Allison  (1978) 
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Table  17.  Relative  Photoionization  Intensities  for  Acetylene 
Excited  by  Mg  K  (1253.6  eV). 


Orbital 

F0a 

(2) 
S(E) 

(2) 
S^E^SHIFT 

17.79 

0PWb 
26.86 

r       b 
Exp. 

2a 

g 

18.90 

17.75 

13.3 

3°9 

0.51 

0.48 

0.48 

0.65 

not  observed 

l\ 

1.00 

1.00 

1.00 

1.00 

1.00 

2a 
u 

7.74 

- 

7.50 

11.18 

9.5 

Frozen  orbital  approximation. 
bCavell  and  Allison  (1978). 
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Table  18.  Relative  Photoionization  Intensities  for  Acetylene 
Excited  by  He(II)  (40.81  eV). 


Orbital 

F0a 
0.45 

(2) 
Z(E) 

0.36 

(2) 
E^SHIFT 

0.72 

0PWb 
0.09 

Exp. 

2°g 

0.26 

3°g 

0.45 

0.44 

0.87 

0.77 

0.68 

Uu 

1.00 

1.00 

1.00 

1.00 

1.00 

2c 
u 

0.29 

- 

0.50 

0.41 

0.29 

Frozen  orbital  approximation. 
bCavell  and  Allison  (1978). 
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The  relative  intensities  in  Table  17  computed  with  the  diagonal 
2p-h  TDA  and  Mg  K  source  have  been  plotted  in  Fig.  6.  As  in  Fig.  5, 
the  lines  have  been  given  Lorentzian  widths  and  the  experimental  spec- 
trum appears  as  an  insert.  The  most  striking  difference  between  these 
two  spectra  is  the  absence  of  the  intense  shake-up  in  the  theoretical 
spectrum.  Not  only  was  this  peak  predicted  to  lie  7  eV  higher  than  the 
experimental  peak  in  energy,  it  also  yielded  a  relative  intensity  several 
orders  of  magnitude  smaller  than  the  2a  intensity  (at  23.4  eV).  The 
weak  experimental  peak  at  about  14.5  eV  arises  from  the  ionization  of 
2a   electrons  by  Mg  K^     .  radiation  (Cavell  and  Allison,  1978)  and  is 
absent  in  the  theoretical  spectrum.  The  3a  peak  which  should  occur  at 
16.8  eV  is  not  observed  experimentally  but  can  be  identified  as  a  shoulder 
on  the  2a  peak  in  the  theoretical  spectrum. 

Figure  7  shows  the  dependence  of  the  photoionization  cross-section 
on  photon  energy  from  0-200  eV.  Over  this  energy  range,  ionization  from 
the  Itt  orbital  is  predicted  to  be  most  probable,  and  this  is  verified 
by  the  experimental  data  of  Cavell  and  Allison  (1978)  obtained  at  21.2  eV 
(He  (  I  )),  40.8  eV  (He  (II)),  and  151.4  eV  (Zr  M  ).  Another  interesting 
feature  of  this  figure  is  the  slight  minimum  exhibited  by  the  2a  curve 
around  125  eV.  Minima  of  this  type  in  the  photoionization  cross-section 
were  first  predicted  by  Cooper  (1962)  and  are  referred  to  as  Cooper 
minima.  Cooper  minima  have  been  observed  in  the  cross-sections  for  a 
number  of  atoms  by  means  of  X-ray  and  ultraviolet  absorption  spectroscopy 
(see  e.g.  Codling,  1976);  however,  at  present  there  are  no  photoelectron 
spectroscopic  data  of  this  type  owing  to  the  limitation  of  photon  sources 
which  are  available.  The  application  of  synchrotron  radiation  to  photo- 
electron  spectroscopy  should  soon  offer  a  means  of  studying  the  energy 


Figure  6.  A  plot  of  the  theoretical  ESCA  spectrum  for  the  valence 
ionizations  of  acetylene.  The  experimental  spectrum  of 
Cavell  and  Allison  (1978)  is  sketched  in  the  insert. 
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Figure  7.  A  plot  of  the  photoionization  cross-sections  versus 
photon  energy  for  the  valence  orbitals  of  acetylene 
in  the  region  0-200  eV. 
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dependence  of  the  orbital  cross-sections  in  a  continuous  energy  interval 
from  0  to  about  200  eV  (Codling,  1973). 

Orbital  plots  are  presented  for  the  2a   ,  3a  ,  and  Itt  Feynman-Dyson 

u    g       u        J 

amplitudes  in  Fig.  8.  Since  the  dominant  component  in  each  of  these 
amplitudes  is  the  corresponding  Hartree-Fock  orbital  (i.e.  2o  ,  3o  ,  and 
Itt  ,  respectively)  the  correlation  and  relaxation  corrections  are  not 
readily  observable.  In  order  to  examine  these  many-electron  contributions 
more  readily,  density  difference  plots  between  the  amplitudes  of  Fig.  8 
and  their  corresponding  Hartree-Fock  orbital  were  made  and  are  presented 
in  Figs.  9-11.  Since  the  Feynman-Dyson  amplitudes  are  not  normalized, 
it  is  necessary  to  normalize  them  before  computing  the  density  difference. 
In  all  of  these  plots,  the  Hartree-Fock  density  is  negative  which  means 
any  positive  distortions  imply  density  enhancement  in  the  propagator  am- 
plitude while  negative  distortions  imply  density  diminution.   Figure  9 
shows  a  density  diminution  in  both  the  C-C  and  C-H  bonding  regions  with 
a  density  enhancement  in  the  anti-bonding  region.  Figure  10  shows  an 
enhancement  in  the  C-C  pi  bonding,  and  Fig.  11  shows  an  enhancement  in 
the  C-H  bonding  with  a  slight  diminution  of  anti-bonding  character. 

It  is  apparent  from  the  numerical  results  presented  in  this  section 
that  the  calculated  relative  photoionization  intensities,  at  least  within 
the  plane  wave  approximation,  are  not  very   sensitive  to  improvements  in 
the  Feynman-Dyson  amplitudes.   It  is  likely  that  more  accurate  photoelec- 
tron  wavefunctions  will  improve  the  sensitivity  of  these  quantities  but 
not  at  the  orthogonalized  plane  wave  level.  Although  the  OPW  approxima- 
tion is  more  justifiable  formally,  the  results  of  Rabalais  et  aj_.  (1974) 
and  Cavell  and  Allison  (1978)  do  not  exhibit  any  marked  improvements  to 
the  plane  wave  results  in  the  two  molecules  studied  here.  Qualitative 


Figure  8.     Orbital    plots   for  the  2c?u,   3ag,   and   1ttu   Feynman-Dyson 
amplitudes  of  acetylene.     Approximate  scales  on  the 
plots  are  ±0.72  a.u.   for  2ou,  ±0.28  a.u.   for  3aa,  and 
±0.29  a.u.    for   1ttu.  y 
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Figure  9.  A  density  difference  plot  between  the  3aq  Feynman-Dyson 
amplitude  and  the  3og  Hartree-Fock  orbital  of  acetylene. 
Each  contour  represents  a  density  increment  of  about 
1.5  x  10-4  atomic  units. 
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Figure  10.  A  density  difference  plot  between  the  1ttu  Feynman-Dyson 
amplitude  and  the  1ttu  Hartree-Fock  orbital  of  acetylene. 
Each  contour  represents  a  density  increment  of  about 
3.6  x  10"3  atomic  units. 
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Figure  11.  A  density  difference  plot  between  the  2au  Feynman-Dyson 
amplitude  and  the  2au  Hartree-Fock  orbital  of  acetylene. 
Each  contour  represents  a  density  increment  of  about 
4.4  x  10_4  atomic  units. 
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changes  in  the  propagator  amplitudes  may  be  observed  via  density 
difference  plots;  however,  more  quantitative  comparisons  must  await 
other  computational  applications. 


CONCLUSIONS  AND  EXTENSIONS 

The  primary  objective  of  this  dissertation  has  been  the  investiga- 
tion of  alternative  decouplings  which  may  allow  more  accurate  calcula- 
tions of  the  electron  propagator  with  less  computational  expenditure 
than  the  operator  product  decouplings.  Contrary  to  previous  contentions 
that  decouplings  of  the  propagator  equations  of  motion  and  choice  of 
reference  state  averages  are  interrelated  (Oddershede  and  J0rgensen, 
1977),  we  have  demonstrated  that  they  are  in  fact  independent  approxima- 
tions, hence  justifying  the  systematic  improvement  of  one  or  the  other. 
We  have  chosen  to  examine  the  decoupling  approximation  while  maintaining 
a  Hartree-Fock  reference  state  since  these  reference  states  have  a  simple 
theoretical  representation  and  are  efficiently  generated  by  a  number  of 
available  computer  programs.  Since  Hartree-Fock  reference  states  ignore 
electron  correlation  (by  definition),  the  correlation  and  relaxation 
corrections  to  the  electron  propagator  must  be  exclusively  incorporated 
through  the  decoupling  approximation. 

The  superoperator  formalism  which  was  originally  introduced  as  a 
notational  convenience  has  a  rich  algebraic  structure  that  has  been 
largely  unexploited.  As  we  demonstrate  in  Chapter  3,  the  superoperator 
Hamiltonian  may  be  partitioned  and  a  perturbation  theory  may  be  developed 
in  the  operator  space.  The  perturbation  expansion  of  the  superoperator 
resolvent  readily  yields  the  Dyson  equation  and  allows  an  identification 
of  wave  and  reaction  superoperators .  The  definition  of  these  auxiliary 
superoperators  elucidates  the  parallelism  between  the  superoperator 
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formalism  and  the  diagrammatic  expansion  method  and  permits  a  unified 
conceptual  framework. 

When  the  superoperator  resolvent  is  approximated  by  an  inner  pro- 
jection, various  choices  of  the  inner  projection  manifold  which  corre- 
spond to  different  decouplings  may  be  viewed  as  different  approximations 
to  the  wave  superoperator  in  the  following  equation  (Born  and  Ohrn,  1977] 

|h)  =  W(E)  |a)  . 
In  particular,  the  operator  product  decoupling  corresponds  to  the  choice 

W  =  I  +  {a^} 

~   a    +       -j-  + 
W  =  I  +  {a,  a,}  +  {a,1  a, a  a  } 
k  1     u  ran 

the  moment  conserving  decoupling  corresponds  to 

W  =  I  +  H 

W  =  I  +  H  +  H2 
*  .  .  , 

and  the  diagram  conserving  decoupling  of  the  Pade'  type  corresponds  to 

W(E)  =  I  +  TQ(E)V 

W(E)  =  I  +  T0(E)VTQ(E)V 

Alternatively,  the  superoperator  resolvent  may  be  approximated  by 
an  outer  projection  in  which  case  decoupling  approximations  correspond 
to  various  approximations  to  the  reaction  (or  self-energy)  superoperator. 
Two  decouplings  of  this  type  which  have  been  presented  are  the  simple 
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diagram  conserving  decoupling  which  corresponds  to 

t(E)  =  V  +  VT(E)V 
with 

T(E)  =  TQ(E) 

T(E)  =  TQ(E)  +  T0(E)VTQ(E) 

*    5 

and  the  renormalized  decoupling  which  corresponds  to 

t(E)  =  V  +  VT(E)V 
with 

T(E)  =  TQ(E)  +  TQ(E)MT(E) 

and  where  M  is  an  effective  potential  obtained  by  projecting  the  pertur- 
bation superoperator  onto  a  subspace  of  the  full  operator  space.  When 
this  subspace  is  chosen  to  consist  of  simple  operator  products,  this 
decoupling  is  formally  equivalent  to  the  operator  product  decoupling, 
however,  other  approximations  may  be  envisioned. 

Each  of  the  above-mentioned  decoupling  approximations  has  been 
studied  in  detail  computationally  and  has  been  evaluated  on  the  basis 
of  a  comparison  of  the  propagator  poles  with  experimental  ionization 
energies.  The  moment  conserving  decouplings  discussed  in  Chapter  2  were 
found  to  yield  unacceptably  poor  numerical  results.  It  is  plausible 
that  since  the  moment  matrices  involve  averages  of  various  powers  of 
the  full  superoperator  Hamiltonian,  the  moment  expansion  may  not  be  con- 
vergent anywhere  in  the  complex  plane.  If  this  is  the  case,  analytic 
continuation  cannot  be  performed,  and  the  Pade'  approximant  method  will 
not  improve  convergence.  In  any  event,  it  may  be  generally  concluded 
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that  the  number  of  conserved  moments  is  no  indication  of  the  accuracy 
of  this  decoupling. 

A  more  successful  approach  is  the  diagram  conserving  decouplings 
presented  in  Chapter  3.  A  simple  second-order  self-energy  truncation 
was  found  to  overcorrect  the  Hartree-Fock  orbital  energy  estimate  of 
the  ionization  energy  and  yielded  results  that  were  generally  too  small. 
The  inclusion  of  third-order  diagrams,  particularly  rings  and  ladders,  is 
necessary  to  obtain  reasonable  agreement  with  experiment.  Energy  shifts 
resulting  from  the  addition  of  diffuse,  polarization  functions  to  the 
computational  basis  were  observed  to  be  nearly  the  same  magnitude  as  the 
shifts  obtained  from  third-order  ring  and  ladder  diagrams,  hence  care 
must  be  taken  to  insure  basis  saturation. 

The  renormalized  2p-h  TDA  decoupling  was  found  to  be  the  most  satis- 
factory decoupling  studied.  Although  only  a  diagonal  approximation  was 
implemented,  principal  ionization  energies  accurate  to  within  an  electron 
volt  were  generally  obtained.  Shake-up  energies  were  less  accurately 
described  owing  to  the  neglect  of  ion  state  relaxation  and  hole-hole 
correlation  in  the  2p-h  TDA. 

The  numerical  results  of  Chapter  5  indicate  that  the  calculation  of 
relative  photoionization  intensities  using  a  plane  wave  approximation 
for  the  photoelectron  is  not  a  sensitive  reflection  of  the  quality  of 
the  Feynman-Dyson  amplitudes.  The  severity  of  the  plane  wave  approxima- 
tion most  likely  obliterates  the  many-electron  relaxation  and  correlation 
corrections.  More  accurate  photoelectron  representations  such  as  orthog- 
onal ized  plane  waves  or  Coulomb  waves  should  afford  a  more  accurate  eval- 
uation of  the  Feynman-Dyson  amplitudes,  or  alternatively,  one-electron 
properties  may  be  computed  using  the  single-particle,  reduced  density 
matrix  obtained  from  the  spectral  density. 
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There  are   several  other  possible  extensions  of  this  research  that 
should  be  mentioned.  Most  obvious  is  an  application  of  these  decoupling 
ideas  to  the  polarization  propagator.  Currently,  decouplings  of  this 
propagator  are  based  on  operator  products.  The  odd  or  Fermion-1 ike 
operator  products  in  the  electron  propagator  theory  must  be  replaced  by 
even  or  Boson-like  operator  products  to  account  for  the  particle-conserv- 
ing nature  of  the  polarization  propagator  theory,  but  with  this  minor 
modification,  the  application  of  the  diagram  conserving  and  renormalized 
decouplings  also  seems  feasible. 

More  accurate  decouplings  are  still  needed  for  the  electron  propa- 
gator in  order  to  better  describe  shake-up  energies.  With  this  goal  in 
mind,  Cederbaum  and  co-workers  have  recently  implemented  the  full,  non- 
diagonal  2p-h  TDA  by  exploiting  spin  and  orbital  symmetry  to  reduce  the 
dimension  of  the  2p-h  operator  subspace.  Preliminary  results  indicate 
that  this  approximation  provides  some  improvement  but  is  still  inadequate. 
The  need  for  higher  operator  products,  particularly  quintuple  products, 
has  been  demonstrated  by  Bagus  and  Viinikka  (1977)  in  CI  studies  and  may 
necessitate  the  implementation  of  a  3p-2h  renormal ization. 

Another  attractive  avenue  for  extension  is  a  combination  of  the 
2p-h  TDA  renormalization  with  the  Hamiltonian  partitioning  of  Kurtz  and 
Ohrn.  Some  work  along  this  line  is  in  progress  by  Ortiz  and  Ohrn  (1979) 
and  preliminary  results  are  yery   encouraging.  One  criticism  of  this 
approximation  is  that  if  one  interprets  it  as  an  alternative  partitioning 
of  the  superoperator  Hamiltonian  in  which  the  eigenvalues  of  the  unper- 
turbed Hamiltonian  correspond  to  AE(SCF)  energies  (as  in  Eq.  (3.23)), 
then  the  denominators  of  the  self-energy  which  are  obtained  from 

( EI-HJL )   (some  operator  product) 
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should  also  involve  AE(SCF)  energies.  For  example,  the  second-order 
self-energy  denominator  would  be 


[E  +  AEk(SCF)  -  AE^SCF)  -  AE^SCF)]"1 


Since  it  is  probably  not  feasible  to  perform  AE(SCF)  calculations  for 
all  doublet  (N-l)  and  (N+l)-electron  states,  simple  approximations  to 
the  AE(SCF)  energies  such  as 


AEk(SCF)   s  E|< 


E 
a,p 


<akj  lpkxpk|  |ak> 


e  -  e 
a    p 


or 


<ak 


AE.(SCF)  ^  e.  -  I     7 r— 

k        k   a.p^^p^^ 


pk><pk| |ak> 


pa>  +  <pk|  |pk>  -  <akMak> 


might  be  considered. 

Finally,  the  problem  of  open  shell  reference  states  has  not  been 
considered  here.  The  retention  of  a  restricted  Hartree-Fock  reference 
state  in  this  regard  will  be  a  tremendous  simplification,  and  the  various 
spin  and  orbital  symmetry  couplings  between  the  reference  state  and 
operator  manifold  may  be  treated  with  the  theory  of  tensor  operators. 
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